Source: https://pubs.rsc.org/en/content/articlehtml/2019/cp/c8cp04445g?page=search
Timestamp: 2019-04-26 06:21:31+00:00

Document:
Although the signatures of natural optical activity have been recognized for over two centuries1 their fruitful application for the direct determination of absolute stereochemical configuration has been revolutionized only relatively recently by the advent of reliable quantum-chemical paradigms for predicting such chiroptical properties from first principles.2 Of particular interest is the dispersive phenomenon of circular birefringence (CB),3,4 whereby the differential retardation of helicity components experienced by an electromagnetic wave as it traverses an isotropic ensemble of chiral molecules allows for the facile discrimination of enantiomers through the opposing directions (yet equal angles) of induced polarization rotation. Linear-response calculations built upon coupled cluster (CC) and density functional theory (DFT) frameworks have proven remarkably successful at reproducing the intrinsic behavior observed for rigid species under rarefied (vapor-phase) conditions,5,6 often yielding quantitative agreement with experiments when purely electronic contributions dominate the specific optical rotation (SOR), [α]Tλ (in deg dm−1 (g mL−1)−1), measured at a given non-resonant wavelength (λ) and near-ambient temperature (T). For more complex systems entrained in chemically relevant condensed phases (e.g., solutions), the accord between experiment and theory tends to deteriorate markedly owing to a conspiracy of competing factors that can alter both the magnitude and the sign optical rotatory dispersion (ORD or wavelength-resolved CB),7 including environmental perturbations (solute–solvent and/or solute–solute coupling), conformational flexibility (large-amplitude nuclear displacements), and vibrational motion (small-amplitude nuclear dynamics). The incorporation of such effects into a comprehensive theoretical description of optical activity imparts formidable challenges and represents an active area of research. The present work strives to explore the boundaries of current computational models by analyzing wavelength-resolved rotatory powers acquired for dilute aqueous solutions of non-rigid (R)-glycidyl methyl ether (R-GME), where complications incurred from each of the aforementioned difficulties are compounded further by the formation of strong hydrogen bonds among adjoining solute and solvent molecules. To be more precise, the focus of this study is not to achieve the best possible agreement between theory and experiment, but to explore how the myriad of factors contributing to the overall ORD of a complex system in aqueous solution are compounded and to demonstrate that, despite the notable success achieved by various computational protocols in the past, substantial challenges remain to be addressed before a reliable approach can be formulated to correctly capture the essential chemistry and physics governing this surprisingly elusive property (viz. producing accurate numerical results for the right reasons).
When considering the quantitative simulation of fundamental molecular properties such as the natural optical activity of a complex system in the condensed phase, many different types of effects should be taken into account and, to this day, a definitive computational protocol for accomplishing such analyses does not exist. The role of the quantum-mechanical model (CC, DFT, or other) and the accompanying basis set must be evaluated carefully alongside more intricate environmental effects arising from inevitable interactions between the target (solute) and its surroundings (solvent). This difficulty is compounded further by the fact that dispersive chiroptical response depends strongly on the details of molecular structure and, therefore, non-rigid species that readily explore vast regions of conformation space can prove particularly difficult to model.
As suggested above, the highly structured conformation landscapes that mediate the spectral signatures evoked from flexible chiral substrates are expected to be affected markedly by strong solute–solvent interactions – a scenario epitomized by the pervasive hydrogen-bonding networks that characterize liquid water. Consequently, the ensuing work has adopted several distinct theoretical approaches to model the dispersive chiroptical properties observed for dilute aqueous solutions of R-GME. Initial conformer-averaging results have been obtained from a refined implementation of the PCM ansatz16–18 that properly includes crucial effects arising from local-field corrections, nonequilibrium solvation, and vibrational motion, the latter being evaluated at anharmonic levels of theory.19 To overcome the fundamental limitations inherent to such continuum-dielectric methods, subsequent computational efforts have exploited a multiscale explicit/implicit solvation formalism whereby reliable quantum-mechanical (QM) calculations of natural optical activity are coupled with accurate molecular-mechanical (MM) treatments of environmental perturbations. These analyses were based on the recently developed mixed quantum-classical QM/FQ scheme,20,21 which accounts for important mutual polarization phenomena in the classical MM portion (viz., the cybotactic region proximate to the solute) through use of fluctuating charges (FQ). More specifically, by averaging the QM/FQ rotatory powers predicted independently for a series of solute–solvent configurations emerging from extensive molecular dynamics (MD) simulations, the pronounced influence of conformational and vibrational degrees of freedom can be incorporated automatically without needing to rely on otherwise ad hoc assumptions (e.g., individual conformers residing in putatively isolated potential-energy minima). In all cases, theoretical estimates of wavelength-resolved specific rotation have been compared with their experimental counterparts to elaborate the microscopic origins of R-GME chiroptical properties as measured under solvated conditions. The diverse models deployed during the present investigation have permitted various aspects of solvation to be examined systematically, including factors that alter the structural and response properties of flexible solutes capable of forming hydrogen bonds with the surrounding solvent. Through careful analysis of the effects that different computational protocols can exert on emerging properties, we highlight the fact that resorting to crude approximations can lead to apparent agreement with experiments; however, such observations may be deceptive and often represent a fortuitous cancellation of errors that might not be transferable to other systems. As such, our findings warn against the temptation of applying simple solutions to a problem that remains, to this day, quite complex and in need of further exploration.
Requisite samples of (R)-glycidyl methyl ether (R-GME) were obtained from a commercial source (Sigma-Aldrich) at 97% stated purity and used without further treatment. Although the percentage enantiomeric excess (%ee) of this material was not listed by the manufacturer, values of specific optical rotation (SOR) measured for sodium D-line (589.3 nm) excitation were consistent with the highest magnitudes of this quantity previously published.22 As such, observed rotatory powers were not subjected to post-collection corrections designed to account for chemical and/or enantiomeric imperfections.
Solution-phase studies of dispersive optical activity were performed in an ambient, room-temperature (∼25°) quartz sample cell (10.000 ± 0.005 cm length) through use of a commercial polarimeter (Perkin-Elmer 341; ±0.002° angular accuracy) operating at discrete visible and ultraviolet excitation wavelengths filtered from Na I and Hg I atomic-emission lamps. The solvent employed for the present investigation of aqueous R-GME was of spectrophotometric grade (Alfa-Aesar), with solute concentrations being kept as low as possible (typically <10−3 g mL−1 or <10−2 M) so as to maintain superior angular precision while minimizing deleterious effects arising from possible solute–solute interactions (e.g., aggregation). The latter processes were quantified by performing a series of dilution experiments, the details of which have been relegated to the ESI.† Experimental values of wavelength-resolved rotatory power compiled in Table 1 represent the results of linear least-squares regressions designed to extrapolate solution-phase measurements to the limit of infinite dilution so as to determine the intrinsic (solvated) response of the chiral solute.23 To further elucidate the provenance of dispersive chiroptical signatures, an electronic circular dichroism (ECD) spectrum of aqueous R-GME was acquired by exploiting a commercial instrument (Applied Photophysics Chirascan; 0.2 nm resolution); however, the lowest-lying absorption feature was found to reside beyond the solvent cut-off wavelength of roughly 185 nm, thereby precluding any detailed assessment of pertinent excitation energies.
The absence of atomistic detail (for the solvent) in the QM/PCM treatment precludes the quantitative modeling of specific solute–solvent interactions such as those known to arise in the presence of aqueous hydrogen bonding, and fails to account for dynamical processes taking place between solvent and solute molecules, which can cause the latter to explore regions of conformational space that might otherwise be considered inaccessible. To overcome these limitations, we employed a polarizable QM/MM approach in which each solvent atom is treated classically, but is endowed with a fluctuating charge (FQ) that responds dynamically to the solute electrostatic potential.32–35 Requisite fluctuating-charge values were determined through systematic parameterization by assigning a hardness and an electronegativity to each type of atom.36,37 The coupling between the QM and MM layers introduces explicit solvent-dependent terms into the quantum-mechanical Hamiltonian which then propagate into all response properties and result in explicit solvation effects.21 Analogous QM/FQ schemes have been deployed successfully to evaluate diverse molecular properties and spectroscopic signatures,21,38–47 including those of ORD;20,48 however, the basic FQ methodology also can be used in MD simulations alone.49 Calculations of properties that involve the charge moments beyond the electric dipole, including the magnetic dipole and thus ORD, can suffer from a gauge-origin dependence of the results.50–57 This problem is addressed by means of gauge-including atomic orbitals,50,51 which also can be used in conjunction with the FQ method.41 The force field governing the classical trajectories of aqueous R-GME was created through application of the Joyce program58 by performing a series of reference DFT calculations at the B3LYP/6-311G(d) level of theory. This work utilized the conductor-like PCM (C-PCM)59 to incorporate water interactions implicitly and strove to reproduce key intramolecular features, with particular attention directed towards the distributions of large-amplitude torsional coordinates (vide infra).
Lahiri et al.8 have reported significant differences between SOR values predicted for R-GME by means of DFT and the Linear Response Coupled-Cluster Singles and Doubles (LR-CCSD) scheme, where respective use of gauge-including atomic orbitals (GIAOs)50,51 and the modified velocity-gauge (MVG) ansatz52 ensured translational invariance. Notwithstanding the significant effects that can be incurred by the action of the surrounding medium, a rigorously correct treatment of chiroptical properties should not prescind the intrinsic (electronic) behavior exhibited by the unsolvated molecule. Building upon optimized conformational geometries and relative Boltzmann weights obtained from DFT calculations (vide infra), Fig. 1 contrasts ORD curves computed for isolated (vapor-phase) R-GME at the B3LYP and LR-CCSD levels of theory by employing the aug-cc-pVDZ and aug-cc-pVTZ basis sets. The resulting profiles are quite similar at longer wavelengths, but deviate significantly for λ < 350 nm. The origins of natural optical activity can be envisioned through an infinite sum-over-states expression where each successive term entails the rotatory strength of an electronically excited manifold divided by a factor that scales in proportion to the detuning of the incident optical frequency from that of the attendant electronic transition. As such, the relative ORD contributions made by the lowest-lying electronic state can be expected to grow rapidly as the applied wavelength is tuned into resonance. From this perspective, the disparity among the curves in Fig. 1 can be considered to reflect differences in estimated first excitation energies, which will affect the spectral position of the associated pole in the response function for the system. This assertion is confirmed in Fig. 2 where a direct comparison of energies calculated for the lowest-lying electronic state of each R-GME structural isomer reveals LR-CCSD estimates uniformly to be larger than their B3LYP counterparts by roughly 10% (corresponding to 0.6 eV or 20 nm). The histograms in Fig. 1 also show only modest differences to exist between B3LYP/aug-cc-pVTZ and B3LYP/aug-cc-pVDZ results, reinforcing that a correlation-consistent basis of double-ζ quality affords sufficient computational reliability. Similar analyses could be performed by incorporating the PCM description of solvation; however, in contrast to the DFT QM/MM treatments that form the crux of the present study, a comparable implementation of the LR-CCSD/FQ methodology does not exist. Consequently, ensuing efforts will concentrate on the putative roles of solvation and on the various DFT-based theoretical approaches that can be deployed to elucidate their ramifications for SOR. This will abrogate attempts to achieve a full quantitative agreement with experiment, which otherwise would demand the combination of more expensive electronic-structure methods along with the sophisticated solvation models elaborated below.
Fig. 1 Calculated optical activity in vacuo. B3LYP/aug-cc-pVDZ, B3LYP/aug-cc-pVTZ, and LR-CCSD/aug-cc-pVDZ predictions of ORD profiles are depicted for isolated R-GME. These results follow from conformer-averaging procedures based on geometries and relative energies optimized at the B3LYP/aug-cc-pVDZ level of theory. While DFT analyses of optical activity utilized the length gauge with origin-invariance guaranteed by use of GIAOs, their LR-CCSD counterparts employed the MVG ansatz.
Fig. 2 Calculated excitation energies in vacuo. B3LYP/aug-cc-pVDZ, B3LYP/aug-cc-pVTZ, and LR-CCSD/aug-cc-pVDZ predictions of the first vertical excitation energy are shown for each R-GME conformer. These results stem from ground-state geometries optimized at the B3LYP/aug-cc-pVDZ level of theory.
Building upon the isolated-molecule picture, the solvation environment can be introduced systematically though use of the PCM ansatz. This leads to the set of nine structural isomers depicted in Fig. 3, each of which was tested to represent a minimum of the Potential Energy Surface (PES) by the absence of imaginary (harmonic) vibrational frequencies. The dihedral angles (θ and ϕ) and relative free-energies (ΔG) extracted from emerging B3LYP/aug-cc-pVTZ/PCM equilibrium geometries are reported in Table 2 along with their computed Boltzmann population fractions at T = 298 K (f298K). The similarity of these parameters to those previously reported8 for PCM calculations of R-GME in acetonitrile underscores the inability of the implicit-solvation approach to distinguish among highly polar media.
Fig. 3 R-GME conformational flexibility. The nine stable ground-state minima predicted for aqueous R-GME by QM/PCM analyses are illustrated, with each conformation being distinguished by a Roman numeral as well as by approximate values of dihedral angles θ and ϕ. The latter coordinates are defined explicitly on the accompanying geometrical structure, which also establishes the numbering scheme used to identify heavy atoms (adapted from Lahiri et al.8).
Aside from the relative populations of contributing species, successful application of the conformer-averaging ansatz to the dispersive optical activity of aqueous R-GME demands that chiroptical signatures be evaluated for each equilibrium configuration of the solvated molecular framework. The reliability of this treatment was assessed by exploiting a variety of different functionals and the aug-cc-pVDZ basis set to perform DFT linear-response calculations of SOR, with Fig. 4 highlighting the sodium D-line (589.3 nm) rotatory powers obtained by implementing four popular approaches (cf. ESI† for other wavelengths and functionals). The depicted specific-rotation parameters are rather consistent, displaying neither alterations in sign nor significant changes in magnitude for an individual conformer. The robustness of these results affords confidence that our analyses can be extended to the polarizable QM/MM paradigm without worrying that conclusions drawn about solvent-induced effects might be affected strongly by the nature of the underlying DFT methodology. Indeed, given the similarity of SOR predictions obtained from various functionals, ensuing analyses have focused exclusively on B3LYP calculations. At shorter wavelengths, more dramatic disparities among the predictions obtained with different functionals are inevitable since each places the poles for the attendant response function (viz. resonances) at unique locations. Such factors can lead to pronounced changes in ORD (both magnitude and sign) as the wavelength decreases, with ensuing analyses addressing these effects by including results from several functionals.
Fig. 4 Conformer-specific optical activity. Values of sodium D-line (589.3 nm) SOR predicted for each low-lying stable conformation of aqueous R-GME (as optimized at the B3LYP/aug-cc-pVDZ/PCM level of theory) are depicted for four choices of the density functional used in DFT/PCM linear-response calculations.
where λ0 signifies the putative location of a low-lying excited state that contributes in proportion to the amplitude factor A and the final term accounts for the cumulative effects of other electronic resonances far removed from the excitation wavelength (λ0 ≪ λ). Detailed analyses of the solution-phase experimental data yield A = 1.20(12) × 107 deg dm−1 (g mL−1)−1 nm2, B = −1.92(13) × 107 deg dm−1 (g mL−1)−1 nm2, and λ0 = 174.7(46) nm (one standard deviation uncertainties on final significant digits in parentheses), thereby reinforcing the assertion that the lowest-lying electronic manifold of aqueous R-GME resides below the solvent cutoff wavelength of ∼185 nm. Analogous least-squares regressions performed on linear-response computations allowed emerging SOR parameters to be interpolated between and extrapolated beyond discrete wavelength values.
Fig. 5 Ensemble averaged optical activity. The experimental ORD profile for aqueous R-GME is compared with conformer-averaging predictions made at the B3LYP/aug-cc-pVDZ/PCM level of theory. In addition to purely electronic calculations (based on stable minimum-energy geometries), the results obtained by applying anharmonic vibrational corrections at temperatures of 0 K and 298 K are shown.
When comparing calculations to measurements, one should keep in mind that the stated goal of the present work is to highlight the possible effects that various aspects of the computational protocol can exert on evaluated properties, as opposed to seeking a unique combination of theoretical methods (perhaps solely applicable to the R-GME system) that maximizes the level of agreement attained. Nevertheless, a direct juxtaposition of theoretical and experimental findings does serve as an indispensable guide for evaluating both the magnitude and the direction of trends emerging from such benchmarking efforts.
The overall shape of the observed ORD curve for aqueous R-GME is well reproduced by QM/PCM calculations, although the calculated SOR rises more rapidly than experiment as the incident optical frequency increases (cf.Fig. 5). This behavior might reflect the inability of our quantum-chemical method to reproduce the position of the first pole in the DFT response function, which (when underestimated) can cause the resulting dispersion profile to shift upwards at shorter wavelengths. Such assertions are supported by Fig. 2, where vapor-phase predictions of excitation energies by LR-CCSD methods were found to be uniformly larger than their DFT counterparts. Unfortunately, the intrinsic absorption cutoff of water hampered efforts to measure the location of the lowest-lying electronic state for this system, thus precluding a direct comparison between measured and calculated excitation energies. One also finds that vibrational corrections seem to enhance the accord between theory and experiment; however, this improvement is not sufficient to predict the correct sign of dispersive optical activity at all wavelengths. Another source of discrepancy might reside in the relative populations of individual conformers. As is evident from Fig. 4, the various structural isomers of solvated R-GME possess antagonistic chiroptical properties that can differ in both magnitude and sign. Consequently, even small inaccuracies in conformational sampling could lead to a rather large variations in ensemble-averaged properties. Finally, the solvation model itself could be an important source of error, since the use of a continuum ansatz perforce neglects all effects arising from the atomistic nature of the solvent, including the hydrogen-bonding motifs that are expected to play a crucial role in aqueous media. To address the latter two points, subsequent analyses have enlisted the more refined polarizable QM/MM approach to capture the essential details of solute–solvent coupling.
2.2.3 Explicit solvation: polarizable QM/MM based on FQs. As made clear by the continuum-solvation analyses of the previous section, the most challenging aspect of simulating the dispersive optical activity for aqueous R-GME resides in the concurrent existence of nine low-lying structural isomers, each of which possesses vastly different chiroptical properties and fractional populations that, in turn, are influence strongly by environmental perturbations. This problem is exacerbated by the fact that rotation about the dihedral angles θ and ϕ (cf.Fig. 3), which gives rise to the distinct geometries considered by the canonical conformer-averaging ansatz, is not impeded sufficiently to consider the attendant local minima of the PES as representing truly isolated entities in solution. Consequently, the entire conformational landscape spanned by these large-amplitude nuclear degrees of freedom should be sampled to give a more realistic description of the solvated system – a goal that can be achieved readily though use of classical molecular dynamics (MD) simulations.
To enhance theoretical treatments of coupling between R-GME and the surrounding solvent, two additional interaction sites (each designated as a virtual site or VS) have been added to the molecular topology of the solute. In particular, these entities have been inserted at the centroid positions of the molecular orbitals related to the lone pairs of the oxygen atoms, the essential characteristics of which were determined through use of Boys method60 with the resulting geometrical features allowing for the permanent electric-dipole moment of R-GME to remain essentially unchanged. Our implementation of virtual sites yields a much more chemically accurate description of hydrogen bonding among solute and solvent molecules by enabling the correct directionality for key O⋯H–O binding motifs to be incorporated from the onset. Indeed, similar approaches have been shown to improve the results of MD simulations systematically in the case of hydrogen-bonding media,61,62 although the use of such virtual sites within a polarizable QM/MM theoretical framework remains largely unexplored. It should be remarked that virtual sites are not of a fundamental nature, and that the solute–solvent interactions responsible for the directionality of hydrogen bonds can indeed originate from underlying quantum-mechanical interactions, which cannot be captured fully using a classical model. In this regard, virtual sites constitute an attempt to sample the correct solute–solvent conformational space without resorting to computationally prohibitive quantum dynamics simulations.
In order to accurately parametrize the force field (FF) needed to describe large-amplitude torsional degrees of freedom, a set of 48 and 24 relaxed R-GME structures were computed about the θ and ϕ dihedral angles, respectively. The results of this procedure are summarized in Fig. 6, where good agreement is realized between the emerging FF predictions and the DFT reference calculations. A cut of the PES along the θ coordinate exhibits a global minimum-energy configuration at θ ≈ −90°, which is separated by a substantial barrier of nearly 25 kJ mol−1 from another metastable geometry at θ ≈ 45°. A third low-lying conformation is found in the vicinity of θ ≈ 180°. The analogous section of the PES taken along ϕ is quasi-symmetric about the point ϕ = 0 with the minimum-energy structure at ϕ ≈ 180°, which maximizes displacement of the side-chain methyl group from the epoxide ring, being isolated by barriers of nearly 10 kJ mol−1 from higher-energy forms located at ϕ ≈ ±80° (cf.Fig. 3).
Fig. 6 Force field description. The molecular force field (FF; solid red curve) for aqueous R-GME, as deduced by implementing a fitting procedure with the Joyce program, is compared with the results of relaxed quantum-mechanical calculations (QM; discrete black symbols) performed at the B3LYP/6-311G(d)/C-PCM level of theory. The two panels highlight the torsional potentials obtain by separately varying the large-amplitude dihedral angles θ (on the left) and ϕ (on the right).
A 50 ns molecular-dynamics trajectory (NVT ensemble) was performed for a system composed of a lone R-GME solute and roughly 2800 water molecules, the latter being described by the TIP3P model.63 A total of 8000 snapshots were extracted from the last 40 ns of the MD simulation (one snapshot every 5 ps), with initial analyses of these results focusing on the specificity of intermolecular interactions. As shown in Fig. 7, quantitative information regarding the nature of hydrogen-bond (HB) formation can deduced from the radial distribution functions, g(r), evaluated between the oxygen atoms of water and the O1(epoxide)/O2(side-chain) atoms of R-GME (cf.Fig. 3 for numbering scheme).
Fig. 7 Vestiges of solvation-shell structure. Radial distribution functions, g(r), resulting from MD trajectories of aqueous R-GME that include virtual sites (MD2) are depicted as a function of the distance, r, between oxygen atoms. Two computed profiles are shown, highlighting displacements between either O1 (continuous blue curve) or O2 (dashed green curve) of the chiral solute molecule and the oxygen center of the water solvent.
Even a cursory inspection of the depicted g(r) profiles reveals several key differences that reflect the subtle nature of solute–solvent coupling. The O1–OH2O distribution shows an intense peak of 1.25 height at 0.26 nm, clearly suggesting the presence of strong hydrogen bonding. The coordination number of 2.03 computed for the first minimum of this g(r) curve indicates that the heterocyclic oxygen couples, on average, with two water molecules. Solvent-mediated interactions appear to be weaker in the case of O2, where a coordination number of 1.4 has been determined for the first hydration shell. Such disparities might stem from differential accessibility of the two R-GME oxygen centers. Whereas O1, with virtual sites jutting out from the epoxide ring, can be approached readily by water, both the mobility of the side-chain and the bulk of the nearby methylene/methyl groups may introduce steric hindrance about the analogous acceptor sites of O2. The second hydration shell for O2 appears to be somewhat more pronounced than that of O1, as documented by the corresponding g(r) amplitudes of 1.32 and 1.21, respectively.
The hydrogen-bonding network of aqueous R-GME was examined further by counting the effective number of non-covalent linkages as a function of time. Here, a simple geometric criterion was adopted whereby each OR-GME–OH2O interaction was defined as being a hydrogen bond if (i) the computed distance between the two atoms was less than 0.35 nm, and (ii) the computed ∠HH2OOH2OOR-GME angle was less than or equal to 30°. Each frame extracted from the final 5 ns of the MD trajectory has been considered for this analysis.
The resulting histograms for the computed number of hydrogen bonds are illustrated in Fig. 8 and essentially confirm the insights deduced from the g(r) radial-distribution profiles. Most of the inspected configurations reveal two non-covalent linkages adjoining O1, while nearly 70% of them display a solitary hydrogen bond involving O2. It is worth noticing that a non-negligible number of frames suggest three-fold coupling to O1, thereby corroborating the predominant role that this epoxide oxygen center plays for interactions with the solvent residues.
Fig. 8 Extent of hydrogen bonding. The effective number of hydrogen bonds formed to the O1 (left panel) and O2 (right panel) oxygen atoms of the aqueous R-GME solute are depicted in the form of histograms, as computed from the final 5 ns of the MD simulation including virtual sites (MD2).
A more comprehensive analysis of hydrogen-bonding motifs can be found in the spatial distribution functions (SDFs) reported in Fig. 9, which highlight findings obtained from tandem MD simulations performed without (MD1) and with (MD2) the presence of virtual sites (VSs) on the solute oxygen atoms. Such three-dimensional maps show the regions of space where solvent molecules congregate, thereby affording insights into the differences that exist among intermolecular binding patterns. Given that R-GME explores numerous configurations during an MD trajectory, these figures were generated by considering the solvent positions about O1 or O2 relative to the plane formed by this reference oxygen atom and the two adjoining (bonded) carbons, with the emerging SDF being superimposed on a representative solute geometry. While significant disparities do not exist between the MD1 and MD2 treatments for O2, marked variations can be seen in the analogous O1 results. More specifically, the absence of VS interactions (MD1) for O1 tends to produce crown-like layers of high probability above this lone center of the epoxide ring. The incorporation of O1 virtual sites (MD2) alters this scenario to yield two distinct areas of solvent occupation, the disposition of which reflects the sp3 hybridization of the solute oxygen atom and the expected directionality of attendant hydrogen bonds. Such contrasting behavior can be expected to influence computed chiroptical properties (vide infra).
Fig. 9 Spatial distribution functions. Probability density maps are shown for the water solvent molecules closest to the oxygen atoms of the R-GME solute. These results stem from the MD1 (top row) and MD2 (bottom row) simulations, with the probability distributions for aqueous oxygen and hydrogen centers being color coded in red and white, respectively.
Theoretical analyses of dispersive optical activity for aqueous R-GME utilized a series of 8000 solute–solvent configurations emerging from the aforementioned MD trajectories. For each resulting “snapshot,” a sphere of 12 Å radius fixed about the geometrical center of the solute was extracted and the wavelength-resolved SOR was calculated through use of the polarizable QM/FQ scheme by exploiting the SPC FQ parameters given by Rick et al.32 Simulations were performed both without (MD1) and with (MD2) virtual sites on the solute oxygen atoms and the configuration-averaged chiroptical response at each discrete wavelength was evaluated through use of four different DFT functionals. Fig. 10 presents a comparison between the predicted ORD curves for aqueous R-GME and their experimental counterpart.
Fig. 10 QM/FQ simulations of optical activity. The experimental ORD profile for aqueous R-GME is contrasted with QM/FQ predictions of chiroptical response as obtained through the use of four different density functionals for linear-response calculations. The two panels illustrate the influence of excluding (top) and including (bottom) virtual sites on the oxygen atoms of the chiral solute.
For both sets of MD simulations (MD1 and MD2) shown Fig. 10, the placement of ORD curves generated by the various DFT functionals are in the same relative order, such that B3LYP and M06-2X respectively predict the highest and lowest rotatory powers at any given wavelength. More importantly, significant disparities clearly exist between the MD1 and MD2 results, with the former uniformly underestimating experimental measurements while the latter consistently overestimates them. The incongruity between these models can be attributed to two primary effects, one pertaining to the solute and the other related to the attendant solvation shell. The presence of virtual sites on oxygen centers clearly influences the types of interactions that can take place between the solute and solvent molecules. This can produce different conformational distributions for R-GME that, as demonstrated by the QM/PCM findings of the previous section, can alter the magnitude and/or sign of the ensemble-averaged chiroptical response markedly. However, even if relative conformer populations were unaffected, the arrangement of solvent molecules around the solute can be modified significantly by virtual sites, thus leading to pronounced changes in observed properties. The veracity of the latter claim is supported by the SDF maps of Fig. 9, which highlight the differential localization of solvent molecules obtained from the MD1 and MD2 treatments.
Fig. 11 Configurational dependence of optical activity. The sodium D-line (589.3 nm) rotatory powers predicted for individual solute–solvent configurations emerging from MD simulations of aqueous R-GME are depicted by discrete markers that are color coded according to the optical rotation scale in the legend. These results follow from QM/FQ calculations implemented through use of the B3LYP density functional. The two panels illustrate the influence of excluding (top) and including (bottom) virtual sites on the oxygen atoms of the chiral solute, leading to average SOR parameters of −32 and −7 deg dm−1 (g mL−1)−1, respectively, as compared to the experimental value of −17 deg dm−1 (g mL−1)−1.
Fig. 12 Conformational landscape for optical activity. Predictions of sodium D-line (589.3 nm) rotatory powers obtained at the B3LYP/aug-ccpVDZ/PCM level of theory are depicted for individual low-lying conformers of aqueous R-GME. For each structural isomer, a circle color-coded according to optical-rotation value (see legend) has been centered at the optimized values of θ and ϕ dihedral coordinates with its area being proportional to the computed QM/PCM Boltzmann population fraction at T = 298 K. The corresponding conformer-averaged SOR parameter of −24 deg dm−1 (g mL−1)−1 must be compared to the experimentally measured result of −17 deg dm−1 (g mL−1)−1.
The results emerging from our MD simulations of aqueous R-GME unequivocally demonstrate that dispersive optical activity is inherently a dynamical quantity. Especially in cases of conspicuous conformational flexibility, theoretical treatments based on distinct energy minima do not afford a realistic physical description of the solvated species, even when corrections for small-amplitude vibrational displacement have been included. Furthermore, the strong dependence of SOR on the precise nature of intermolecular interactions suggests that an accurate modeling of the composite solute–solvent system is essential (e.g., even for configurations supporting nearly identical geometries of the solute molecule), with conventional continuum models potentially failing to provide the level of detail needed to reliably characterize this property.
While one might argue that the ORD curves in Fig. 10b are quite similar to those in Fig. 5, such observations cannot be taken as evidence that discrete and continuum models are interchangeable. Indeed, dynamical processes enable the solute molecule to explore a much broader structural landscape than that associated with isolated isomeric forms, with each region of conformational space contributing uniquely to the overall (ensemble-averaged) chiroptical response in a manner that can deviate significantly from corresponding QM/PCM results. As has been shown for very similar systems,20 QM/FQ and QM/PCM calculations may yield dramatically different predictions. Thus, it appears to be fortuitous (viz. an exception rather than a rule) that the combined changes in conformational dynamics and electronic-response properties afforded by the QM/FQ method produces only small alterations in computed ORD profiles for aqueous R-GME relative to the QM/PCM scheme.
This study has exploited theoretical calculations and experimental measurements to unravel the complex relationships that exist between the dispersive optical activity of a conformationally flexible chiral molecule, (R)-glycidyl methyl ether (R-GME), and the pervasive interactions sustained from a surrounding environment of liquid water. For the targeted system, our efforts have revealed that the traditional conformer-averaging ansatz, which relies on the presumption of independent and isolated structural isomers “solvated” by means of implicit continuum models, is inadequate for the quantitative description of circular-birefringence (CB) phenomena. A refined approach to treat solute–solvent coupling, the mixed quantum-classical polarizable QM/MM scheme based on fluctuating charges (FQs),20,21 was enlisted to predict the wavelength dependence of specific optical rotation (SOR) for aqueous R-GME and to elucidate the provenance of observed chiroptical signals.
Of paramount importance to the present work has been the structural and electronic effects incurred by an aqueous medium, with particular emphasis directed towards the marked influence that microsolvation can exert on conformational dynamics and response properties. The emerging results highlight how seemingly minor details of the computational protocol, such as the presence or absence of virtual sites on solute oxygen atoms, can have enormous repercussions for predicted chiroptical behavior, despite the fact that the utilized QM/FQ approach still considers the encompassing solvent as a perturbation acting on the chiral solute (viz., only the latter is treated fully quantum mechanically and thus is responsible for generating the molecular response). This also helps to explain why continuum treatments may be inappropriate for describing the SOR of aqueous systems that can form strong hydrogen bonds. In particular, the PCM and kindred models are unable to account for the directionality of intermolecular interactions – a crucial feature that is recovered in the atomistic QM/FQ paradigm. However, even in the case of MD simulations, the intrinsic specificity of hydrogen-bond linkages cannot be realized without the incorporation of virtual sites that have been shown to dramatically alter conformational distributions and chiroptical signatures. Indeed, the pronounced variations in the magnitude and sign of computed rotatory powers for R-GME as the configuration of the surrounding water evolves in space and time highlights how observed macroscopic phenomena arise from a complex amalgamation of contributions made by vastly different microscopic species.
The present results show that much work remains to be done before an accurate and transferable computational model for the ORD of flexible molecules in aqueous solution can be realized. While approximate methods like those relying on continuum models of solvation may give predictions in reasonable accord with experimental measurements for some systems, this observation might stem from a fortuitous cancellation of errors arising from an unbalanced description for the solute conformational landscape and the hydrogen-bonding patterns formed with surrounding solvent. Particular care must be taken in such cases, with specific analyses being performed to examine critical aspects of solute–solvent coupling and to verify the key assumptions underlying more simplistic treatments.
Finally, it must be emphasized that the conclusions reached from the present investigations of optical activity do not necessarily carry over to the modeling of other quantities. Dispersive chiroptical signatures appear to be exceptionally sensitive to the instantaneous configurations of both solute and solvent molecules, as reflected by their strong dependence on subtle characteristics of the microsolvation environment. Because SOR stems from a highly non-localized electronic response property that is not easily separated into contributions arising from different portions of the system, its accurate prediction for solvated species is especially challenging and definitive answers regarding the physical origins of such phenomena have proven elusive. Further work clearly is needed to identify and characterize the dominant paradigms that govern solvation processes in diverse surroundings, including non-polar media for which the inclusion of electrostatic interactions alone may not be sufficient. Recent progress in the realistic theoretical description of such non-electrostatic effects has been reported for the energetics of solute–solvent coupling,64 though formal extension to the realm of molecular-response functions still is under development.
We are thankful for the computer resources provided by the high performance computer facilities of the SMART Laboratory (http://smart.sns.it/). CC gratefully acknowledges the support of H2020-MSCA-ITN-2017 European Training Network “Computational Spectroscopy In Natural sciences and Engineering” (COSINE), grant number 765739. We also would like to thank Matteo Ambrosetti for his help in generating the images shown in Fig. 9. P. M. L. and P. H. V. gratefully acknowledge the support of the U.S. National Science Foundation under the auspices of grant CHE-1464957.
T. M. Lowry, Optical Rotatory Power, Longmans, Green and Company, London, 1935 Search PubMed.
Comprehensive Chiroptical Spectroscopy, ed. N. Berova, P. L. Polavarapu, K. Nakanishi and R. W. Woody, John Wiley & Sons, Hoboken, NJ, 2012 Search PubMed.
L. Barron, Molecular Light Scattering and Optical Activity, Cambridge University Press, Cambridge, 2nd edn, 2004 Search PubMed.
P. H. Vaccaro, Comprehensive Chiroptical Spectroscopy, John Wiley & Sons, Hoboken, NJ, 2012, ch. 11, pp. 275–323 Search PubMed.
P. Lahiri, K. B. Wiberg, P. H. Vaccaro, M. Caricato and D. T. Crawford, Angew. Chem., Int. Ed., 2014, 53, 1386–1389 CrossRef CAS.
P. Lahiri, K. B. Wiberg and P. H. Vaccaro, J. Phys. Chem. A, 2017, 121, 8251–8266 CrossRef CAS.
T. Aharon, P. Lemler, P. H. Vaccaro and M. Caricato, Chirality, 2017, 30, 383–395 CrossRef.
P. Lahiri, K. B. Wiberg and P. H. Vaccaro, J. Phys. Chem. A, 2015, 119, 8311–8327 CrossRef CAS PubMed.
T. Müller, K. B. Wiberg and P. H. Vaccaro, J. Phys. Chem. A, 2000, 104, 5959–5968 CrossRef.
T. Müller, K. B. Wiberg, P. H. Vaccaro, J. R. Cheeseman and M. J. Frisch, J. Opt. Soc. Am. B, 2002, 19, 125–141 CrossRef.
S. M. Wilson, K. B. Wiberg, J. R. Cheeseman, M. J. Frisch and P. H. Vaccaro, J. Phys. Chem. A, 2005, 109, 11752–11764 CrossRef CAS.
Continuum Solvation Models in Chemical Physics, ed. B. Mennucci and R. Cammi, John Wiley & Sons Ltd, Chichester, England, 2007 Search PubMed.
J. Tomasi, B. Mennucci and R. Cammi, Chem. Rev., 2005, 105, 2999–3093 CrossRef CAS.
B. Mennucci, WIREs Comput. Mol. Sci., 2012, 2, 386–404 CrossRef CAS.
P. Mukhopadhyay, G. Zuber, P. Wipf and D. N. Beratan, Angew. Chem., Int. Ed., 2007, 46, 6450–6452 CrossRef CAS PubMed.
C. Cappelli, F. Lipparini, J. Bloino and V. Barone, J. Chem. Phys., 2011, 135, 104505 CrossRef.
B. Mennucci, C. Cappelli, R. Cammi and J. Tomasi, Chirality, 2011, 23, 717–729 CrossRef CAS.
S. Pipolo, R. Cammi, A. Rizzo, C. Cappelli, B. Mennucci and J. Tomasi, Int. J. Quantum Chem., 2011, 111, 826–838 CrossRef CAS.
F. Egidi, J. Bloino, V. Barone and C. Cappelli, J. Chem. Theory Comput., 2012, 8, 585–597 CrossRef CAS.
F. Lipparini, F. Egidi, C. Cappelli and V. Barone, J. Chem. Theory Comput., 2013, 9, 1880–1884 CrossRef CAS.
C. Cappelli, Int. J. Quantum Chem., 2016, 116, 1532–1542 CrossRef CAS.
W. Felzmann, D. Castagnolo, D. Rosenbeiger and J. Mulzer, J. Org. Chem., 2007, 72, 2182–2186 CrossRef CAS.
P. L. Polavarapu, A. Petrovic and F. Wang, Chirality, 2003, 15, S143–S149 CrossRef CAS.
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 16 Revision A.03, Gaussian Inc., Wallingford CT, 2016 Search PubMed.
T. Yanai, D. P. Tew and N. C. Handy, Chem. Phys. Lett., 2004, 393, 51–57 CrossRef CAS.
Y. Zhao and D. G. Truhlar, Theor. Chem. Acc., 2008, 393, 215–241 Search PubMed.
C. Adamo and V. Barone, J. Chem. Phys., 1998, 110, 664–675 CrossRef.
S. W. Rick, S. J. Stuart and B. J. Berne, J. Chem. Phys., 1994, 101, 6141–6156 CrossRef CAS.
S. W. Rick and B. J. Berne, J. Am. Chem. Soc., 1996, 118, 672–679 CrossRef CAS.
S. W. Rick, S. J. S. J. S. Bader and B. J. Berne, J. Mol. Liq., 1995, 65–66, 31–40 CrossRef.
F. Lipparini and V. Barone, J. Chem. Theory Comput., 2011, 7, 3711–3724 CrossRef CAS PubMed.
W. J. Mortier, K. van Genechten and J. Gasteiger, J. Am. Chem. Soc., 1985, 107, 829–835 CrossRef CAS.
R. Chelli and P. Procacci, J. Chem. Phys., 2002, 117, 9175–9189 CrossRef CAS.
F. Egidi, R. Russo, I. Carnimeo, A. D'Urso, G. Mancini and C. Cappelli, J. Phys. Chem. A, 2015, 119, 5396–5404 CrossRef CAS.
T. Giovannini, M. Olszówka and C. Cappelli, J. Chem. Theory Comput., 2016, 12, 5483–5492 CrossRef CAS.
T. Giovannini, M. Olszówka, F. Egidi, J. R. Cheeseman, G. Scalmani and C. Cappelli, J. Chem. Theory Comput., 2017, 13, 4421–4435 CrossRef CAS PubMed.
F. Lipparini, C. Cappelli and V. Barone, J. Chem. Phys., 2013, 138, 234108 CrossRef PubMed.
I. Carnimeo, C. Cappelli and V. Barone, J. Comput. Chem., 2015, 36, 2271–2290 CrossRef CAS PubMed.
F. Lipparini, C. Cappelli and V. Barone, J. Chem. Theory Comput., 2012, 8, 4153–4165 CrossRef CAS PubMed.
F. Lipparini, C. Cappelli, G. Scalmani, N. D. Mitri and V. Barone, J. Chem. Theory Comput., 2012, 8, 4270–4278 CrossRef CAS.
M. Caricato, F. Lipparini, G. Scalmani, C. Cappelli and V. Barone, J. Chem. Theory Comput., 2013, 9, 3035–3042 CrossRef CAS.
T. Giovannini, G. Del Frate, P. Lafiosca and C. Cappelli, Phys. Chem. Chem. Phys., 2018, 20, 9181–9197 RSC.
F. Egidi, G. Lo Gerfo, M. Macchiagodena and C. Cappelli, Theor. Chem. Acc., 2018, 137, 82 Search PubMed.
F. Egidi, I. Carnimeo and C. Cappelli, Opt. Mater. Express, 2015, 5, 196–209 Search PubMed.
G. Mancini, G. Brancato and V. Barone, J. Chem. Theory Comput., 2014, 10, 1150–1163 CrossRef CAS PubMed.
F. London, J. Phys. Radium, 1937, 8, 397–409 CrossRef CAS.
J. R. Cheeseman, M. J. Frisch, F. J. Devlin and P. J. Stephens, J. Phys. Chem. A, 2000, 104, 1039–1046 CrossRef CAS.
T. B. Pedersen, H. Koch, L. Boman and A. M. J. S. de Méras, Chem. Phys. Lett., 2004, 393, 319–326 CrossRef CAS.
L. N. Vidal, F. Egidi, V. Barone and C. Cappelli, J. Chem. Phys., 2015, 142, 174101 CrossRef.
L. N. Vidal, T. Giovannini and C. Cappelli, J. Phys. Chem. Lett., 2016, 7, 3585–3590 CrossRef CAS.
K. Ruud, T. Helgaker and P. Bouř, J. Phys. Chem. A, 2002, 106, 7448–7455 CrossRef CAS.
P. J. Lestrange, F. Egidi and X. Li, J. Chem. Phys., 2015, 143, 234103 CrossRef.
S. Bernadotte, A. J. Atkins and C. R. Jacob, J. Chem. Phys., 2012, 137, 204106 CrossRef.
V. Barone, I. Cacelli, N. De Mitri, D. Licari, S. Monti and G. Prampolini, Phys. Chem. Chem. Phys., 2013, 15, 3736–3751 RSC.
V. Barone and M. Cossi, J. Phys. Chem. A, 1998, 102, 1995–2001 CrossRef CAS.
S. F. Boys, Rev. Mod. Phys., 1960, 32, 296–299 CrossRef CAS.
M. Macchiagodena, G. Mancini, M. Pagliai and V. Barone, Phys. Chem. Chem. Phys., 2016, 18, 25342–25354 RSC.
M. Macchiagodena, G. Mancini, M. Pagliai, G. Del Frate and V. Barone, Chem. Phys. Lett., 2017, 677, 120126 CrossRef.
P. Mark and L. Nilsson, J. Phys. Chem. A, 2001, 105, 9954–9960 CrossRef CAS.
T. Giovannini, P. Lafiosca and C. Cappelli, J. Chem. Theory Comput., 2017, 13, 4854–4870 CrossRef CAS.

References: V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V. 

V. 

V. 
 V. 
 V.