CPL Spectra of Camphor Derivatives in Solution by an Integrated QM/MD Approach

1. Introduction

Circularly polarized luminescence (CPL) and electronic circular dichroism (ECD) are complementary techniques giving access to the properties of the excited and ground electronic states, respectively for chiral chromophores or chromophores embedded in chiral environments (Longhi et al., 2016; Tanaka et al., 2018). However, the magnetic dipole transition moment is generally much smaller than the electric one, thus hindering the investigation of systems showing strong luminescence (Riehl and Muller, 2012). Therefore, for both chiroptical spectroscopies, electric dipole forbidden and magnetic dipole allowed transitions represent a particularly appealing situation.

Camphor and its derivatives, and in particular their formally forbidden n → π^* carbonyl transitions, have been extensively studied by chiroptical methods, from both experimental and computational points of views (Dekkers and Closs, 1976; Schippers et al., 1983; Pritchard and Autschbach, 2010; Longhi et al., 2013; McAlexander and Crawford, 2015; Duong and Fujiki, 2017). Nevertheless, simulation of the camphor CPL is far from being straightforward because of (i) the inherent dependence of the vibrational information from solvation interactions and (ii) the bisignate nature of the spectra. This latter phenomenon was recently qualitatively investigated by Longhi et al. (2013) and Duong and Fujiki (2017). According to Kasha's rule (Kasha, 1950), luminescence is expected to occur from the lowest excited electronic state. Therefore, within the Franck-Condon approximation, the first ECD and CPL transitions should have the same sign. In fact, the bisignate shape is a consequence of the pyramidalization of the carbonyl moiety, which gives rise to two low-energy conformers connected through a large amplitude motion (LAM) in the first excited state. This interconversion path needs to be properly accounted in order to reproduce and evaluate the different vibronic contributions in the spectra (Cerezo et al., 2016; Baiardi et al., 2017).

In this contribution, we apply a QM/MD approach to perform accurate simulations of camphor one-photon absorption (OPA) and emission (OPE) spectra and of their chiral counterparts ECD and CPL. We employ our recently developed ONIOM/EE-PMM (Del Galdo et al., 2019, 2020) procedure to model solvent effects on the electronic properties of the chromophore. The method belongs to the general category of multiscale procedures and merges variational and perturbative approaches in the calculation of spectroscopic features. Starting from extensive classical simulations of a target electronic state of a chromophore, effective clustering procedures are applied in order to identify a relevant set of sub-trajectories within the complete sampling. Within each cluster, one reference snapshot is utilized for an ONIOM calculation performed within the Electronic Embedding framework [ONIOM/EE, i.e., the variational approach (Dapprich et al., 1999; Vreven et al., 2006; Chung et al., 2015)]. Then, the local fluctuations within each sub-trajectory are treated via the Perturbed Matrix Method [PMM, i.e., the perturbative approach (Aschi et al., 2001; Del Galdo et al., 2018; Zanetti-Polzi et al., 2018)]. The procedure permits to strongly reduce the computational costs when compared to more conventional QM/MM approaches, without lowering the accuracy of the final results.

On top of this, the vibrational contributions to the final spectra are also modeled. We exploit for the purpose a general computational tool that allows the simulation of different kinds of one-photon spectroscopies and supports inclusion of mode-mixings, as well as Franck-Condon (FC) and Herzberg-Teller (HT) effects. Finally, internal coordinates are used in order to minimize the coupling between normal modes and LAMs, which are properly accounted for and treated separately.

2. Methods

2.1. Combining QM and MM Methods, Step 1: Sampling Solute Internal Motions and Solvent Effects

Molecular Dynamics simulations of the ground and first excited electronic states of camphor in methanol were employed to generate the ensembles needed to evaluate the solvatochromic shifts on absorption and emission spectra. In view of the structural rigidity characterizing the camphor ground state, the corresponding MD simulation was run by constraining the solute at its equilibrium geometry. In this way, we did not take into account internal motions through classical approaches, thus adding their contribution a posteriori by means of full QM calculations of vibrational modulation effects by means of models based on the Franck-Condon principle (see sections 2.3). Conversely, as already reported (Longhi et al., 2013; Duong and Fujiki, 2017), in its first electronic state camphor undergoes an internal deformation associated to the carbonyl out of plane bending. After the characterization of the LAM (vide infra), we accounted for its effects by extracting the most representative conformers characterizing the excited state potential energy surface as a function of the carbonyl out of plane bending (namely the two minima of the curve reported in Figure 1). Then, both conformers were utilized for the classical MD sampling performed by constraining the solute geometry while the solvent molecules followed their unconstrained motion. The other internal degrees of freedom of the solute were taken into account in a separate step by full QM computations of vibronic effects within the harmonic approximation. For each MD ensemble the corresponding emission spectra were computed, thus obtaining the final spectra by proper weighting of the results. In Table 1, the relative stabilities of the two conformers are reported, which were used to compute the respective Boltzmann population at 300 K. Remembering that the solute was kept rigid during the MD runs for both ground and excited electronic states, the only properties differentiating the two simulation ensembles were the geometry of the solute and the atomic charges, the latter possibly including virtual sites (VSs) to represent oxygen lone pairs. The positions of VSs were determined by first locating the O-VS distance and the C-O-VS angle for the centroids of localized molecular orbitals at the sp2 oxygen atom of formaldehyde and propanone by using the Boys localization procedure, as reported in Macchiagodena et al. (2016) and Del Galdo et al. (2020). By averaging the (nearly identical) results, we obtained a value of 0.3067 Å for the distance and 111.3 degree for the angle. For the ground electronic state we imposed that the VSs, the O, the C and the other two atoms bound to C (i.e., the two C atoms) lie on the same plane (O-C(O)-C plane). Instead, for the excited state the VSs were constrained in the plane orthogonal to the plane defined by the normal to the C-C-C plane and oxygen atom. Separation of the charge between the oxygen nucleus and the VSs (having vanishing masses) was achieved by imposing that both VSs have the same charge and reproducing the classical dipole moment of the molecule in the target electronic state. The role of VSs was evaluated by performing two MD simulations (with and without their inclusion) for the ground electronic state. Specific solute-solvent interactions were analyzed by computing the probability of hydrogen bonding (HB) between the carbonyl oxygen and the hydroxyl group of methanol with the help of the so-called F function (Pagliai et al., 2003; Del Galdo et al., 2019, 2020), based on exponential decays for deviations from both radial and angular optimal values obtained from pure methanol simulations (we obtained 1.8 Å and 6 degrees for the optimum distance and angle, respectively and 0.2 Å and 5 degrees for the corresponding half width at half maximum, HWHM).

FIGURE 1

Figure 1. Superimposed structures of the two most stable conformers (Conf.1 in blue and Conf.2 in red) and of the transition state between them in the first excite state.

TABLE 1

Table 1. Relative Electronic Energy (ΔE), Free Energy at 300 K (ΔG), Free Energy at 300 K in Solution (ΔG PCM) Values (kJ mol⁻¹) for the first excited state. Electronic energies computed at Cis(d) B2 levels of theory, vibrational and PCM correction computed at TD-DFT B3 level of theory^a.

2.2. Combining QM and MM Methods, Step 2: The ONIOM/EE-PMM Method

Solvatochromic shifts on the different spectra of camphor were evaluated by the ONIOM/EE-PMM procedure we recently proposed (Del Galdo et al., 2019, 2020). The procedure merges variational and perturbative approaches to get a cheap yet accurate characterization of the perturbing effects exerted by an embedding environment on the quantum mechanical properties of a chromophore. The first step is the identification of a set of basins/clusters within a production MD run. For each sub-sampling, a reference frame is chosen as the most representative one of the corresponding cluster and its spectrum is evaluated by a variational procedure (ONIOM/EE) in which the QM model system (here the Camphor molecule) is embedded in a set of point charges representing the environment (here the methanol molecules) (Vreven et al., 2006). For all the other frames of each cluster, the fluctuations of environmental effects with respect to the reference configuration are computed a posteriori by a perturbative approach (PMM). In details, for each frame of each sub-trajectory, the system Hamiltonian is written as the diagonal matrix of the eigenvalues of the reference configuration plus a perturbation matrix representing the difference of the electrostatic potential between the considered frame and the reference value. Diagonalization of this matrix provides a set of eigenvalues (electronic states) representing the instantaneous effects of the embedding environment as provided by the MD trajectories. The operator representing the variation between the perturbing effects exerted by the environment in each frame and the reference, is modeled by exploiting the latest development of the PMM procedure, that is by expanding the perturbing electrostatic potential within the atomic region around each atomic center (atom-based expansion) (Zanetti-Polzi et al., 2018). Interested readers can find details of the procedure in Aschi et al. (2001), Zanetti-Polzi et al. (2018), Del Galdo et al. (2019), and Del Galdo et al. (2020) while specific computational details are reported in the following subsections.

2.3. Combining ONIOM/EE-PMM and Vibronic Calculation

The fist step is the simulation of the vibrational modulation (vibronic) effects in electronic spectra using models based on the Franck-Condon principle and a continuum description (PCM) of bulk solvent effects, as described in detail in Barone et al. (2009), Bloino et al. (2010), and Bloino et al. (2016). In particular, we employed the so-called time-independent (TI) approach, in which the band-shape is obtained as the sum of the individual transitions between the vibrational states of the initial and final electronic states. In order to best decouple the effect of the LAM, we employed the Vertical Gradient (VG) model (Baiardi et al., 2016), which requires the vibrational frequencies of the initial state, but only the forces of the final state at the equilibrium geometry of the initial one (Barone et al., 2009; Bloino et al., 2010, 2016).

Once computed the vibronic spectra for the reference structure of each conformer, they need to be combined with the outcome of the ONIOM/EE-PMM procedure to obtain the final spectra.

At the FC level, a simple shift by the transition energy and scaling by the transition moments for each snapshot leads, after proper averaging, to the final spectrum. Furthermore, on the grounds of previous experience (Del Galdo et al., 2020), effects related to environment fluctuations were considered negligible for magnetic transition dipole moments.

Inclusion of HT terms, accounting for the dependence of the transition dipole moments on the nuclei position, would require the computation of the derivatives of the perturbed transition dipole moments for each MD snapshot. Since this model would become too computationally demanding, we assumed, as in previous works (Del Galdo et al., 2020), that the difference between FC and FCHT spectra is explicitly evaluated only for the reference configuration, whereas the effect of solvent fluctuations is taken as half the value explicitly computed for the FC part as described above. In the following sections, the FCHT acronym will be used to indicate the inclusion of both FC and HT terms in the calculations.

2.4. Computational Details

All geometry optimizations and frequency calculations were carried out with a development version of the Gaussian suite of quantum chemical programs (Frisch et al., 2019), whereas all the MD simulations were performed with the GROMACS software (Berendsen et al., 1995). All the MD production runs of ground and first excited electronic states were performed by constraining the chromophore in the corresponding equilibrium geometry as obtained in the framework of the Density Functional Theory (DFT) (Lee et al., 1988) and its Time Dependent extension (TD-DFT) with the hybrid B3LYP functional (Becke, 1988; Lee et al., 1988) with additional empirical dispersion contributions (D3BJ) (Grimme et al., 2011) and the jul-cc-pVDZ (Dunning, 1989; Papajak et al., 2011) basis set (B3). For each simulation, the solute atomic charges were computed in vacuum at the proper equilibrium geometries by means of the CM5 model (Marenich et al., 2012) at the same level of theory. As mentioned above, two simulations of the ground state were performed including or not virtual sites on the oxygen atom to model directional effects related to sp2 lone pairs. Only the model including virtual sites was used for the simulations of the excited electronic state. All atom Lennard-Jones parameters were taken from the OPLS force field (Jorgensen et al., 1996) MD simulations under periodic boundary conditions were performed in the isothermal-isochoric ensemble (NVT) in order to avoid the additional parameters related to isobaric thermostats required by NPT simulations. The integration step was 2 fs and the temperature was kept constant (300 K) by the velocity-rescaling (Bussi et al., 2007) temperature coupling. Hard geometrical parameters were constrained using the LINCS algorithm (Hess et al., 1997). The particle mesh Ewald method (Darden et al., 1997) was used to compute long range interactions with grid search and cut-off radii of 1.1 nm. A cubic simulation box was utilized and the density was calibrated to obtain in the NVT MD simulations a pressure equal, within the noise, to that provided by a corresponding simulation of the pure solvent. The latter was carried out in the NVT ensemble at T = 300 K and imposing a box density equal to the experimental density of pure methanol at standard conditions (24.58 mol/l, Goodwin, 1987). All the production runs were 10 ns long. We analyzed the solvent effects within each simulation, by computing the probability of the carbonyl oxygen to be engaged in hydrogen bonds with methanol molecules as outlined above (that is, by means of the F function). For each production run, we partitioned the sampling into 4 clusters characterized by two, one (with either VS) or none solute-solvent hydrogen bond, respectively. The partitioning was performed using a threshold of 0.7 for the presence of an effective hydrogen bond. Next, for each cluster, the reference configuration used for the ONIOM/EE calculations is taken as a “collective frame” representative of the average configuration of the molecular environment for the corresponding cluster. To this aim, we extracted 30 snapshots sequentially from the trajectory, from each snapshot we cut a sphere of 30 Å centered around the solute and we assembled them into a collective configuration and we assigned to each environmental atom 1/30 of the actual atomic charge. The first 11 electronic states and the complete matrix of the corresponding dipole moments were computed for the 4 reference configurations exploiting the ONIOM/EE model with the QM part described at the TD-DFT/CAM-B3LYP (Yanai et al., 2004) level in conjunction with the jul-cc-pVDZ (Dunning, 1989; Papajak et al., 2011) basis set. For each electronic state, the corresponding atomic charges were also computed at the same computational level using the CM5 recipe. The above results were then utilized to apply the PMM approach for evaluating fluctuation effects within each cluster. Free energies were determined by adding to electronic energies zero-point energy, thermal, and PCM solvation contributions evaluated in the framework of the rigid rotor/harmonic oscillator approximation (Bloino et al., 2012; Cappelli et al., 2012; Mennucci, 2012). Transition energies have been corrected by the differences between TDDFT/CAMB3LYP/jul-cc-pVDZ and Cis(d) B2PLYPD3/jun-cc-pVTZ (Grimme and Neese, 2007; Jacquemin et al., 2009; Ottochian et al., 2020) (B2) results for isolated camphor. The LAM associated with the double-well potential has been characterized in the framework of the ICPH-model (Baiardi et al., 2017), using 251 DVR basis functions and a step size of 4°. Vibronic spectra were simulated starting from the structure and force constants computed taking into account solvent effects by PCM (Cappelli et al., 2012; Mennucci, 2012). OPE and CPL spectra were computed within the so-called time-independent (TI) approach, employing VG|FC and and VG|FCHT models. The default GAUSSIAN parameters of the class-based integral prescreening scheme were employed (${\text{C}}_1^{max}$ = 50, ${\text{C}}_2^{max}$= 30, ${\text{N}}_I^{max}$ = 10⁹). Normal modes connected with the LAM or methyl rotations were removed. Gaussian distribution functions with HWHM of 900, 550, and 450 cm⁻¹ were used as broadening functions for pure electronic ONIOM/EE-PMM, vibronic PCM, and combined ONIOM/EE-PMM vibronic approaches, respectively. Spatial Distribution Functions plotted as isosurfaces and virtual sites were visualized employing the CAFFEINE software (Salvadori et al., 2016; Lazzari et al., 2020).

3. Results

3.1. Characterization of the LAM

As already mentioned in the introduction, the lowest band in the camphor absorption spectrum is associated with the formally forbidden n → π^* carbonyl transition. However, experiments reveal the presence of a bisignate band in the CPL spectrum that can be connected with a breakdown of the Frank-Condon regime or with the involvement of more than one electronic state in the emission process. Since the second transition lies at a significantly higher energy than the first one, this behavior can be attributed, as previously suggested, to the flexibility of the system. Indeed, in its first electronic state camphor undergoes an internal deformation associated to the carbonyl out of plane bending: a LAM connects two minima through a flat transition state (see Figure 1). These kinds of LAMs are generally ill described by harmonic models (Baiardi et al., 2015, 2016), but use of internal coordinates strongly reduces the coupling of the LAM from the other modes, leading, to a good approximation, to an effective one-dimensional problem (Cerezo et al., 2016; Baiardi et al., 2017). The PES associated with the one-dimensional LAM was explicitly computed through a scan along the CO out-of-plane (OOP) coordinate and it was used to solve numerically the vibrational Schrödiger equation with the help of a quasi-variational approach employing the discrete variable method (Light et al., 1985; Bačić and Light, 1989; Colbert and Miller, 1992; Light and Carrington, 2000; Baiardi et al., 2017). In the upper panels of Figure 2 the PES of the first excited state along the LAM is reported together with the associated vibrational levels obtained with the DVR method, whereas in the bottom panels the ground state PES computed along the same path is reported. Moreover, the vibrational wave functions of those states relevant at 300 K (population above 0.5%) are drawn for both S1 and S0. The different colors highlight the dependence on the LAM of the electric and magnetic transition dipole moment and the angle between them. As expected, in the ground state only one conformer, corresponding to a flat CO structure, is observed along the LAM. On the other hand, the results of the excited state reveal the presence of two conformers corresponding to the out-of-plane bending of the CO, which are separated by a relatively high (18 kJ/mol) transition state (Conf. 1 and Conf. 2). It should also be noted that in S1 the accessible vibrational states at 300 K are well-localized and separated into the two wells. As observed before (Longhi et al., 2013; Duong and Fujiki, 2017), the sign of CPL shows a strong dependence on the LAM. However, as it can be seen in the panels of Figure 2, the most relevant variations are connected with the relative orientation between the two transition dipoles and with the magnitude of the electric transition dipole moment, whereas the corresponding magnetic moment is only negligibly affected by the LAM. Vibronic contributions of the other modes were then described by the Vertical Gradient (VG) model, neglecting the contribution of the LAM already accounted for explicitly. This approach is based on the assumption that the other modes rearrange faster than the LAM (Stendardo et al., 2012; Cerezo et al., 2016; Baiardi et al., 2017).

FIGURE 2

Figure 2. Graphical representation of the PES along the CO out-of-plane bending for the S₁ (upper) and S₀ (lower) electronic states, computed at the B3LYP/Jul-cc-pVDZ (D3BJ) level of theory using the ICPH framework. The vibrational levels and wave functions with contribution above 0.5% at 300K, computed using the variational DVR-based approach, are also reported.

Finally, the strong dependence of the electric transition dipole moments on the nuclear position was accounted for by including Herzberg-Teller (HT) contributions.

3.2. Molecular Dynamics Results

The effectiveness of the interaction between camphor and the embedding solvent was at first assessed by computing the Spatial Distribution Function (SDF) of the methanol molecules around the solute for the ground state simulation. In Figure 3, we compare the results of the SDF calculation performed on simulations involving or not VSs. More precisely, we computed the density distributions of (hydroxyl) hydrogen and oxygen atoms in methanol thus addressing the SDF of the partners of probable HB with camphor oxygen. It is apparent that inclusion of VSs increased the (average) number of methanol molecules around oxygen atom. Moreover, these results show that the arrangement of the solvent is slightly asymmetrical, with a more pronounced effect when the VSs are employed. As a consequence (i) we utilized this simulation as the statistical ensemble to evaluate the solvent effects on the absorption spectra and (ii) also for the production MD runs of the (first) excited electronic state we modeled the camphor including VSs. In Figure 4, the SDFs corresponding to the MD runs of the two conformers simulated to account for the first excited state of camphor are reported. The results suggest that the camphor oxygen (both in the ground and excited states) is engaged in HBs with the solvent molecules. In the following, for the sake of brevity, only the results obtained from the excited state simulations are reported while those corresponding to the ground state simulation can be found in the Supplementary Material (see Figures S1, S2). Nevertheless, it is worth noting that very similar results were found from the classical sampling. We computed the distribution of (i) the distance between each VS and the closest methanol hydroxyl hydrogen and (ii) the angle these two form with the oxygen in methanol. The results obtained from both the excited state simulations are reported in Figure 5. Nearly identical profiles are obtained for the two conformers and, within a given conformer, for each VS: while a clear peak appears in the radial distributions, a less pronounced one characterizes the angular distributions thus suggesting that some geometrical hindrance might weaken HB interactions. Next, we analyzed the HB probability for each excited state simulation by computing the F function described in the section 2. The results reported in Figure 6 show that for both conformers, most of the classical sampling is characterized by the absence of any significant HB. However, we also obtained a non-negligible probability of having one HB for each VS. According to our computations, each methanol molecule was considered to interact with either one VS or the other, thus implying that a probability of around 1 for both the interaction sites corresponds to a total number of 2 HBs for the molecule. On the basis of these results, we decided to partition both the excited state trajectories into a set of clusters, each characterized by a different number of HBs occurring between the camphor and the solvent. As shown in Table 2, the largest cluster accounts for the absence of any significant HB (noHB), then we identified two clusters characterized by one HB depending on the VS involved in the interaction (HB-VS1 and HB-VS2) and the fourth cluster includes the remaining snapshots where at least 2 HBs occur (twoHB). As mentioned above, very similar results were found when the ground state simulation was analyzed, thus allowing us to apply the same partitioning procedure to cluster the ground state trajectory (see Table S1 for the results).

FIGURE 3

Figure 3. Spatial Distribution Function of methanol molecules around camphor from the ground state simulation when the standard atomic charges (left) or the VSs representation (right) are employed. Density distributions of (hydroxyl) hydrogen and oxygen atoms in methanol are represented as white and red volumes, respectively.

FIGURE 4

Figure 4. Spatial Distribution Function of methanol molecules around camphor from the first excited state simulation of Conf. 1 (left) and 2 (right). Density distributions of (hydroxyl) hydrogen and oxygen atoms in methanol are represented as white and red volumes, respectively.

FIGURE 5

Figure 5. Distribution functions of each VS—(methanol hydroxyl) H distance (left) and of each VS—(methanol hydroxyl) H—(methanol hydroxyl) O angles (right) as obtained from the first excited state simulation of camphor in methanol. Upper panels refer to the Conf. 1, lower panels refer to Conf. 2. For all the panels, data computed for the first and second VSs are shown as black and red lines, respectively.

FIGURE 6

Figure 6. Space defined by the values of the F function computed for the each VS in the camphor molecule. Left and right panels correspond to Conf. 1 and Conf. 2, respectively. In each panel the corresponding probability distributions of each F function are also shown.

TABLE 2

Table 2. Population per cent of the clusters within the excited state simulation of Conf. 1 and 2.

The partitioning is the first step toward application of the ONIOM/EE-PMM method, with the four clusters providing the different statistical ensembles from which the representative frames to be employed for the ONIOM/EE calculations were extracted. Then, the PMM procedure was applied to account for fluctuations (within each cluster) of the external electric potential tuning the electronic properties of the chromophore.

3.3. Spectroscopic Results

In order to increase the accuracy of the results, all the transition energies were shifted with reference to the values computed in the gas-phase at the Cis(d) B2PLYPD3/jun-cc-pVTZ level, which is significantly more reliable than TD-DFT with hybrid functionals (Grimme and Neese, 2007; Jacquemin et al., 2009; Laurent and Jacquemin, 2013; Ottochian et al., 2020).

The first step was the simulation of the absorption spectra. The experimental results show an isolated well-defined band for both OPA and ECD spectra without any apparent vibronic sub-structure. According to Longhi et al. (2013), these results suggest a non-negligible interaction with the solvent molecules. In Figure 7, the electronic absorption spectra obtained at the ONIOM/EE-PMM level are reported. Indeed, the four clusters lead to quite different spectra both in terms of energy and intensity, thus contributing to the overall asymmetric shape of both OPA and ECD spectra (see Figure S3 for unweighted cluster spectra). Vibrational modulation effects bring an additional contribution to the final spectra and the corresponding results including bulk solvent effect by mean of the PCM are reported in Figure 8. It is noteworthy that HT terms smooth the prominence of the first band and, especially for the OPA spectrum, cause a non-negligible intensity gain. Finally, as reported in Figure 9, a proper account of specific solvent effects by means of the overall ONIOM/EE-PMM procedure modulates the HT contributions and, in turn, improves the overall agreement with the experimental spectrum, as already observed in the case of the related camphorquinone molecule (Del Galdo et al., 2020).

FIGURE 7

Figure 7. Electronic ECD (left) and OPA (right) spectra of camphor in methanol obtained through the ONIOM/EE-PMM method. The different contributions from each cluster are shown.

FIGURE 8

Figure 8. Vibronic ECD (left) and OPA (right) spectra of camphor in methanol as obtained through VG|FC and VG|FCHT computational procedures in conjunction with PCM model to treat the environment effect.

FIGURE 9

Figure 9. ECD (left) and OPA (right) spectra of camphor in methanol as obtained by combining the ONIOM/EE-PMM procedure outcome with VG|FC (blue line) and VG|FCHT (orange line) models to simulate the vibronic coupling, the corresponding experimental spectrum is also reported (green line). Experimental spectra taken from Longhi et al. (2016).

After having tuned the procedure on absorption spectra we analyzed the corresponding emission spectra. As discussed in section 3.1, the excited electronic state must be described in terms of two disjoint conformers, whose relative stability is a key parameter in the reproduction of the overall spectra. Relying on previous investigations (Barone et al., 2015; Fusè et al., 2019; Paoloni et al., 2020), we employed in this case the B2PLYP-D3 level of theory also accounting for bulk solvent effects by means of the PCM approximation. The results are reported in Table 1, and ΔG PCM values were used to compute Boltzmann population at 300 K in order to get the final simulated spectra.

The electronic emission spectra obtained through application of the ONIOM/EE-PMM method are reported in Figure 10. The differential contributions of both camphor conformers and, for each of them, of the different clusters, are clearly shown. In particular, the presence of different numbers of solute-solvent H-bonds led to quite asymmetrical spectra since these interactions shift the spectrum toward higher energies. On the other hand, the equilibrium between the two main solute conformers is responsible for the inversion of the sign of the transition. In fact, already at the ONIOM/EE-PMM level the bisignate nature of the CPL is well-reproduced. This clearly highlights the relevance of both proper sampling of the solvent distribution and proper weighting of conformer populations. Similarly to the absorption case, the HT term substantially affects the transition intensity. In Figure 11, the OPE and CPL vibronic spectra are reported. It is quite apparent that the CPL intensity of the less populated conformer is higher than that of the prevailing one and it partially compensates the lower population at 300 K. The result of combining vibronic and ONIOM/EE-PMM contributions in the emission spectra is shown in Figure 12. The results are in good agreement with the experimental ones, and it is worth recalling at this point that the full computational procedure does not involve any empirical correction. Moreover, it should be noted that the original CPL spectra from the work of Longhi et al. (2013) are quite noisy (the authors claimed that for both absorption and emission solute-solvent interactions are responsible for the quenching of the vibronic features) and the structure of the experimental CPL signal in Figure 12 is probably due to digitization artifacts. Comparison of the simulated spectra with their experimental counterparts shows that the OPE band-shape is well reproduced by our approach, whereas the agreement is only fair for the CPL spectrum. One of the possible limitations of our computational protocol is the approximate VG vibronic model, which could be too rough to catch and accurately describe the contribution of the vibration on the spectra. Although more refined models could possibly lead to more accurate results, our integrated approach is capable to improve the final results of both ONIOM/EE-PMM and vibronic models taken individually, providing a simulated CPL spectrum significantly closer to the experimental one without any huge increase of the computational cost.

FIGURE 10

Figure 10. Electronic CPL (left) and OPE (right) spectra of camphor in methanol obtained through the ONIOM/EE-PMM method. The different contributions from each cluster for each excited state simulations are shown.

FIGURE 11

Figure 11. Vibronic CPL (left) and OPE (right) spectra of camphor in methanol as obtained through VG|FC and VG|FCHT computational procedures in conjunction with PCM model to treat the environment effect.

FIGURE 12

Figure 12. CPL (left) and OPE (right) spectra of camphor in methanol as obtained by combining the ONIOM/EE-PMM procedure outcome with VG|FC (blue line) and VG|FCHT (orange line) models to simulate the vibronic coupling, the corresponding experimental spectrum is also reported (green line). Experimental spectra taken from Longhi et al. (2016).

4. Conclusions

The main aim of the present contribution was the development of a cheap yet accurate multi-scale strategy for the study of solvatochromic shifts of medium-size flexible chromophores in condensed phases. Explicit treatment of a large amplitude motion and clustering of different MD simulations for each stationary point along the LAM allow the extension of our previous approach from semi-rigid to flexible molecules. At the same time, fitting of atomic charges for excited electronic states permits their MD simulation and the following computation of emission spectra. As a test case, we have chosen a flexible system characterized by remarkable chiroptical properties (namely, the camphor dye in methanol solution), which requires, inter alia, the challenging reproduction of significantly different ECD and CPL spectra. The computational strategy starts from the collection of transition energies and dipole moments tuned by environmental effects provided by MD simulations. Next, vibrational modulation effects on the different electronic transitions are taken into account by a fully quantum mechanical approach including Franck-Condon and Herzberg-Teller contributions together with bulk solvent effects described by a polarizable continuum model. At variance with previous applications, a large amplitude motion has been explicitly considered by means of a discrete-variational approach. Specific chromophore-environment interactions and their fluctuations are described by a recently proposed hybrid procedure, which combines a variational approach (ONIOM/EE), with the PMM perturbative approach. In this way, it becomes possible to obtain accurate results without the need to perform explicit computations for a statistically significant number of solute-solvent configurations. For each stationary point of this flexible system (i.e., the two energy minima of the S1 electronic state and the transition state ruling their inter-conversion corresponding to the single minimum of the ground electronic state) we partitioned the overall simulations into four subsets (clusters) related to the presence or not of strong solute-solvent hydrogen bonds. Although the flexibility of the system increases the number of explicit QM computations with respect to semi-rigid chromophores, the proposed ONIOM/EE-PMM approach strongly reduces the computational effort with respect to a standard QM/MM computation. The good agreement between all the computational results and their experimental counterparts, gives confidence about the robustness and reliability of the overall procedure and permits the analysis of the role of different contributions (stereo-electronic, dynamic, environmental) in tuning the overall spectra.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

Funding

This work has been supported by the Italian MIUR (PRIN 2017, project Physico-chemical Heuristic Approaches: Nanoscale Theory Of Molecular Spectroscopy, PHANTOMS), prot. 2017A4XRCA.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer LE declared a past co-authorship with one of the authors VB to the handling editor.

Acknowledgments

The SMART@SNS Laboratory (http://smart.sns.it) is acknowledged for providing high- performance computer facilities. Prof. Julien Bloino and Dr. Andrea Salvadori are acknowledged for helpful discussions.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem.2020.00584/full#supplementary-material

References

Aschi, M., Spezia, R., Nola, A. D., and Amadei, A. (2001). A first-principles method to model perturbed electronic wavefunctions: the effect of an external homogeneous electric field. Chem. Phys. Lett. 344, 374–380. doi: 10.1016/S0009-2614(01)00638-8

CrossRef Full Text | Google Scholar

Bačić, Z., and Light, J. C. (1989). Theoretical methods for rovibrational states of floppy molecules. Annu. Rev. Phys. Chem. 40, 469–498. doi: 10.1146/annurev.pc.40.100189.002345

CrossRef Full Text | Google Scholar

Baiardi, A., Bloino, J., and Barone, V. (2015). Accurate simulation of resonance-Raman spectra of flexible molecules: An internal coordinates approach. J. Chem. Theory Comput. 11, 3267–3280. doi: 10.1021/acs.jctc.5b00241

PubMed Abstract | CrossRef Full Text | Google Scholar

Baiardi, A., Bloino, J., and Barone, V. (2016). General formulation of vibronic spectroscopy in internal coordinates. J. Chem. Phys. 144:084114. doi: 10.1063/1.4942165

PubMed Abstract | CrossRef Full Text | Google Scholar

Baiardi, A., Bloino, J., and Barone, V. (2017). Simulation of vibronic spectra of flexible systems: Hybrid DVR-harmonic approaches. J. Chem. Theory Comput. 13, 2804–2822. doi: 10.1021/acs.jctc.7b00236

PubMed Abstract | CrossRef Full Text | Google Scholar

Barone, V., Biczysko, M., Bloino, J., Cimino, P., Penocchio, E., and Puzzarini, C. (2015). CC/DFT route toward accurate structures and spectroscopic features for observed and elusive conformers of flexible molecules: pyruvic acid as a case study. J. Chem. Theory Comput. 11, 4342–4363. doi: 10.1021/acs.jctc.5b00580

PubMed Abstract | CrossRef Full Text | Google Scholar

Barone, V., Bloino, J., Biczysko, M., and Santoro, F. (2009). Fully integrated approach to compute vibrationally resolved optical spectra: from small molecules to macrosystems. J. Chem. Theory Comput. 5, 540–554. doi: 10.1021/ct8004744

PubMed Abstract | CrossRef Full Text | Google Scholar

Becke, A. D. (1988). Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A 38, 3098–3100. doi: 10.1103/PhysRevA.38.3098

PubMed Abstract | CrossRef Full Text | Google Scholar

Berendsen, H., van der Spoel, D., and van Druned, R. (1995). GROMACS: a message-passing parallel molecular dynamics implementation. Comput. Phys. Commun. 91, 43–56. doi: 10.1016/0010-4655(95)00042-E

CrossRef Full Text | Google Scholar

Bloino, J., Baiardi, A., and Biczysko, M. (2016). Aiming at an accurate prediction of vibrational and electronic spectra for medium-to-large molecules: an overview. Int. J. Quant. Chem. 116, 1543–1574. doi: 10.1002/qua.25188

CrossRef Full Text | Google Scholar

Bloino, J., Biczysko, M., and Barone, V. (2012). General perturbative approach for spectroscopy, thermodynamics, and kinetics: methodological background and benchmark studies. J. Chem. Theory Comput. 8, 1015–1036. doi: 10.1021/ct200814m

PubMed Abstract | CrossRef Full Text | Google Scholar

Bloino, J., Biczysko, M., Santoro, F., and Barone, V. (2010). General approach to compute vibrationally resolved one-photon electronic spectra. J. Chem. Theory Comput. 6, 1256–1274. doi: 10.1021/ct9006772

PubMed Abstract | CrossRef Full Text | Google Scholar

Bussi, G., Donadio, D., and Parrinello, M. (2007). Canonical sampling through velocity rescaling. J. Chem. Phys. 126:014101. doi: 10.1063/1.2408420

PubMed Abstract | CrossRef Full Text | Google Scholar

Cappelli, C., Bloino, J., Lipparini, F., and Barone, V. (2012). Toward ab initio anharmonic vibrational circular dichroism spectra in the condensed phase. J. Phys. Chem. Lett. 3, 1766–1773. doi: 10.1021/jz3006139

PubMed Abstract | CrossRef Full Text | Google Scholar

Cerezo, J., Mazzeo, G., Longhi, G., Abbate, S., and Santoro, F. (2016). Quantum-classical calculation of vibronic spectra along a reaction path: the case of the ECD of easily interconvertible conformers with opposite chiral responses. J. Phys. Chem. Lett. 7, 4891–4897. doi: 10.1021/acs.jpclett.6b02484

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, L. W., Sameera, W. M. C., Ramozzi, R., Page, A. J., Hatanaka, M., Petrova, G. P., et al. (2015). The ONIOM method and its applications. Chem. Rev. 115, 5678–5796. doi: 10.1021/cr5004419

PubMed Abstract | CrossRef Full Text | Google Scholar

Colbert, D. T., and Miller, W. H. (1992). A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix Kohn method. J. Chem. Phys. 96, 1982–1991. doi: 10.1063/1.462100

CrossRef Full Text | Google Scholar

Dapprich, S., Komaromi, I., Byun, K., Morokuma, K., and Frisch, M. J. (1999). A new ONIOM implementation in Gaussian98. Part I. the calculation of energies, gradients, vibrational frequencies and electric field derivatives. J. Mol. Struct. Theochem. 461–462:1–21. doi: 10.1016/S0166-1280(98)00475-8

CrossRef Full Text | Google Scholar

Darden, T., Tork, D., and Pedersen, L. (1997). Particle mesh Ewald: an n-log(n) method for Ewald sums in large systems. J. Comput. Chem. 18, 1463–1472.

Google Scholar

Dekkers, H. P. J. M., and Closs, L. E. (1976). The optical activity of low-symmetry ketones in absorption and emission. J. Am. Chem. Soc. 98, 2210–2219. doi: 10.1021/ja00424a034

CrossRef Full Text | Google Scholar

Del Galdo, S., Chandramouli, B., Mancini, G., and Barone, V. (2019). Assessment of multi-scale approaches for computing UV-VIS spectra in condensed phases: toward an effective yet reliable integration of variational and perturbative QM/MM approaches. J. Chem. Theory Comput. 15, 3170–3184. doi: 10.1021/acs.jctc.9b00120

PubMed Abstract | CrossRef Full Text | Google Scholar

Del Galdo, S., Fusè, M., and Barone, V. (2020). The ONIOM/PMM model for effective yet accurate simulation of optical and chiroptical spectra in solution: camphorquinone in methanol as a case study. J. Chem. Theory Comput. 16, 3294–3306. doi: 10.1021/acs.jctc.0c00124

PubMed Abstract | CrossRef Full Text | Google Scholar

Del Galdo, S., Mancini, G., Daidone, I., Zanetti Polzi, L., Amadei, A., and Barone, V. (2018). Tyrosine absorption spectroscopy: backbone protonation effects on the side chain electronic properties. J. Comput. Chem. 39, 1747–1756. doi: 10.1002/jcc.25351

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunning, T. H. (1989). Gaussian basis sets for use in correlated molecular calculations. I. the atoms boron through neon and hydrogen. J. Chem. Phys. 90, 1007–1023. doi: 10.1063/1.456153

CrossRef Full Text | Google Scholar

Duong, S. T., and Fujiki, M. (2017). The origin of bisignate circularly polarized luminescence (CPL) spectra from chiral polymer aggregates and molecular camphor: anti-Kasha's rule revealed by CPL excitation (CPLE) spectra. Polym. Chem. 8, 4673–4679. doi: 10.1039/C7PY00958E

CrossRef Full Text | Google Scholar

Frisch, M. J., Trucks, G. W., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., et al. (2019). Gaussian Development Version, Revision J.02. Wallingford, CT: Gaussian, Inc.

Google Scholar

Fusè, M., Mazzeo, G., Longhi, G., Abbate, S., Masi, M., Evidente, A., et al. (2019). Unbiased determination of absolute configurations by vis-à-vis comparison of experimental and simulated spectra: the challenging case of diplopyrone. J. Phys. Chem. B 123, 9230–9237. doi: 10.1021/acs.jpcb.9b08375

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodwin, R. D. (1987). Methanol thermodynamic properties from 176 to 673 k at pressures to 700 bar. J. Phys. Chem. Ref. Data 16, 799–892. doi: 10.1063/1.555786

CrossRef Full Text | Google Scholar

Grimme, S., Ehrlich, S., and Goerigk, L. (2011). Effect of the damping function in dispersion corrected density functional theory. J. Comput. Chem. 32, 1456–1465. doi: 10.1002/jcc.21759

PubMed Abstract | CrossRef Full Text | Google Scholar

Grimme, S., and Neese, F. (2007). Double-hybrid density functional theory for excited electronic states of molecules. J. Chem. Phys. 127:154116. doi: 10.1063/1.2772854

PubMed Abstract | CrossRef Full Text | Google Scholar

Hess, B., Bekker, H., Berendsen, H. J. C., and Fraaije, J. G. E. M. (1997). LINCS: a linear constraint solver for molecular simulations. J. Comput. Chem. 18, 1463–1472. doi: 10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.CO;2-H

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacquemin, D., Wathelet, V., Perpete, E. A., and Adamo, C. (2009). Extensive TD-DFT benchmark: singlet-excited states of organic molecules. J. Chem. Theory Comput. 5, 2420–2435. doi: 10.1021/ct900298e

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Maxwell, D. S., and Tirado-Rives, J. (1996). Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. J. Am. Chem. Soc. 118, 11225–11236. doi: 10.1021/ja9621760

CrossRef Full Text | Google Scholar

Kasha, M. (1950). Characterization of electronic transitions in complex molecules. Discuss. Faraday Soc. 9, 14–19. doi: 10.1039/df9500900014

CrossRef Full Text | Google Scholar

Laurent, A. D., and Jacquemin, D. (2013). TD-DFT benchmarks: a review. Int. J. Quant. Chem. 113, 2019–2039. doi: 10.1002/qua.24438

CrossRef Full Text | Google Scholar

Lazzari, F., Salvadori, A., Mancini, G., and Barone, V. (2020). Molecular perception for visualization and computation: the proxima library. J. Chem. Inf. Model. 60, 2668–2672. doi: 10.1021/acs.jcim.0c00076

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, C., Yang, W., and Parr, R. G. (1988). Development of the colle-salvetti correlation-energy formula into a functional of the electron density. Phys. Rev. B 37, 785–789. doi: 10.1103/PhysRevB.37.785

PubMed Abstract | CrossRef Full Text | Google Scholar

Light, J. C., and Carrington, T. Jr. (2000). “Chapter 4: Discrete-variable representations and their utilization,” in Advances in Chemical Physics, Vol. 114 (New York, NY: John Wiley & Sons, Ltd), 263–310. doi: 10.1002/9780470141731.ch4

CrossRef Full Text | Google Scholar

Light, J. C., Hamilton, I. P., and Lill, J. V. (1985). Generalized discrete variable approximation in quantum mechanics. J. Chem. Phys. 82, 1400–1409. doi: 10.1063/1.448462

CrossRef Full Text | Google Scholar

Longhi, G., Castiglioni, E., Abbate, S., Lebon, F., and Lightner, D. A. (2013). Experimental and calculated CPL spectra and related spectroscopic data of camphor and other simple chiral bicyclic ketones. Chirality 25, 589–599. doi: 10.1002/chir.22176

PubMed Abstract | CrossRef Full Text | Google Scholar

Longhi, G., Castiglioni, E., Koshoubu, J., Mazzeo, G., and Abbate, S. (2016). Circularly polarized luminescence: a review of experimental and theoretical aspects. Chirality 28, 696–707. doi: 10.1002/chir.22647

PubMed Abstract | CrossRef Full Text | Google Scholar

Macchiagodena, M., Mancini, G., Pagliai, M., and Barone, V. (2016). Accurate prediction of bulk properties in hydrogen bonded liquids: amides as case studies. Phys. Chem. Chem. Phys. 18, 25342–25354. doi: 10.1039/C6CP04666E

PubMed Abstract | CrossRef Full Text | Google Scholar

Marenich, A. V., Jerome, S. V., Cramer, C. J., and Truhlar, D. G. (2012). Charge model 5: an extension of hirshfeld population analysis for the accurate description of molecular interactions in gaseous and condensed phases. J. Chem. Theory Comput. 8, 527–541. doi: 10.1021/ct200866d

PubMed Abstract | CrossRef Full Text | Google Scholar

McAlexander, H. R., and Crawford, T. D. (2015). Simulation of circularly polarized luminescence spectra using coupled cluster theory. J. Chem. Phys. 142:154101. doi: 10.1063/1.4917521

PubMed Abstract | CrossRef Full Text | Google Scholar

Mennucci, B. (2012). Polarizable continuum model. WIREs Comput. Mol. Sci. 2, 386–404. doi: 10.1002/wcms.1086

CrossRef Full Text | Google Scholar

Ottochian, A., Morgillo, C., Ciofini, I., Frisch, M. J., Scalmani, G., and Adamo, C. (2020). Double hybrids and time-dependent density functional theory: an implementation and benchmark on charge transfer excited states. J. Comput. Chem. 41, 1242–1251. doi: 10.1002/jcc.26170

PubMed Abstract | CrossRef Full Text | Google Scholar

Pagliai, M., Raugei, S., Cardini, G., and Schettino, V. (2003). Hydrogen bond dynamics in liquid methanol. J. Chem. Phys. 119:6655. doi: 10.1063/1.1605093

CrossRef Full Text | Google Scholar

Paoloni, L., Mazzeo, G., Longhi, G., Abbate, S., Fusè, M., Bloino, J., et al. (2020). Toward fully unsupervised anharmonic computations complementing experiment for robust and reliable assignment and interpretation of IR and VCD spectra from mid-IR to NIR: the case of 2,3-butanediol and trans-1,2-cyclohexanediol. J. Phys. Chem. A 124, 1011–1024. doi: 10.1021/acs.jpca.9b11025

PubMed Abstract | CrossRef Full Text | Google Scholar

Papajak, E., Zheng, J., Xu, X., Leverentz, H. R., and Truhlar, D. G. (2011). Perspectives on basis sets beautiful: Seasonal plantings of diffuse basis functions. J. Chem. Theory Comput. 7, 3027–3034. doi: 10.1021/ct200106a

PubMed Abstract | CrossRef Full Text | Google Scholar

Pritchard, B., and Autschbach, J. (2010). Calculation of the vibrationally resolved, circularly polarized luminescence of D–camphorquinone and (S,S)–trans–β–hydrindanone. ChemPhysChem 11, 2409–2415. doi: 10.1002/cphc.201000054

PubMed Abstract | CrossRef Full Text | Google Scholar

Riehl, J. P., and Muller, G. (2012). “Chapter 3: Circularly polarized luminescence spectroscopy and emission-detected circular dichroism,” in Comprehensive Chiroptical Spectroscopy: Instrumentation, Methodologies, and Theoretical Simulations, Vol. 1 (New York, NY: John Wiley & Sons, Ltd.), 65–90. doi: 10.1002/9781118120187.ch3

CrossRef Full Text | Google Scholar

Salvadori, A., Del Frate, G., Pagliai, M., Mancini, G., and Barone, V. (2016). Immersive virtual reality in computational chemistry: applications to the analysis of QM and MM data. Int. J. Quantum Chem. 116, 1731–1746. doi: 10.1002/qua.25207

PubMed Abstract | CrossRef Full Text | Google Scholar

Schippers, P. H., Van der Ploeg, J. P. M., and Dekkers, H. P. J. M. (1983). Circular polarization in the fluorescence of βγ -enones: distortion in the 1nπ* state. J. Am. Chem. Soc. 105, 84–89. doi: 10.1021/ja00339a015

CrossRef Full Text | Google Scholar

Stendardo, E., Avila Ferrer, F., Santoro, F., and Improta, R. (2012). Vibrationally resolved absorption and emission spectra of dithiophene in the gas phase and in solution by first-principle quantum mechanical calculations. J. Chem. Theory Comput. 8, 4483–4493. doi: 10.1021/ct300664d

PubMed Abstract | CrossRef Full Text | Google Scholar

Tanaka, H., Inoue, Y., and Mori, T. (2018). Circularly polarized luminescence and circular dichroisms in small organic molecules: correlation between excitation and emission dissymmetry factors. ChemPhotoChem 2, 386–402. doi: 10.1002/cptc.201800015

CrossRef Full Text | Google Scholar

Vreven, T., Byun, K. S., Komáromi, I., Dapprich, S., Montgomery, J. A. Jr., Morokuma, K., et al. (2006). Combining quantum mechanics methods with molecular mechanics methods in ONIOM. J. Chem. Theory Comput. 2, 815–826. doi: 10.1021/ct050289g

PubMed Abstract | CrossRef Full Text | Google Scholar

Yanai, T., Tew, D. P., and Handy, N. C. (2004). A new hybrid exchange-correlation functional using the coulomb-attenuating method (CAM-B3LYP). Chem. Phys. Lett. 393, 51–57. doi: 10.1016/j.cplett.2004.06.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Zanetti-Polzi, L., Del Galdo, S., Daidone, I., D'Abramo, M., Barone, V., Aschi, M., et al. (2018). Extending the perturbed matrix method beyond the dipolar approximation: comparison of different levels of theory. Phys. Chem. Chem. Phys. 20, 24369–24378. doi: 10.1039/C8CP04190C

PubMed Abstract | CrossRef Full Text | Google Scholar