Multimodel Approach to the Optical Properties of Molecular Dyes in Solution

TL;DR

This paper presents a multimodel simulation approach combining explicit and implicit solvation models to accurately predict the optical absorption spectra of molecular dyes in solution, validated against experimental data.

Contribution

The study introduces a novel multimodel method that efficiently combines explicit and implicit solvation effects for optical property simulations.

Findings

Time-averaged spectra from multimodel and explicit simulations are highly similar.

The multimodel approach significantly reduces computational resources.

Predicted spectra agree well with experimental measurements.

Abstract

We introduce a multimodel approach to the simulation of the optical properties of molecular dyes in solution, whereby the effects of thermal fluctuations and of dielectric screening on the absorption spectra are accounted for by explicit and implicit solvation models, respectively. Thermal effects are treated by averaging the spectra of molecular configurations generated by an ab initio molecular-dynamics simulation where solvent molecules are treated explicitly. Dielectric effects are then dealt with implicitly by computing the spectra upon removal of the solvent molecules and their replacement with an effective medium, in the spirit of a continuum solvation model. Our multimodel approach is validated by comparing its predictions with those of a fully explicit-solvation simulation for cyanidin-3-glucoside (cyanin) chromophore in water. While multimodel and fully explicit-solvent spectra may differ considerably for individual configurations along the trajectory, their time averages are remarkably similar, thus providing a solid benchmark of the former and allowing us to save considerably on the computer resources needed to predict accurate absorption spectra. The power of the proposed methodology is finally demonstrated by the excellent agreement between its predictions and the absorption spectra of cyanin measured at strong and intermediate acidity conditions.

Full text

\captionsetup

font=small

SISSA]SISSA – Scuola Internazionale Superiore di Studi Avanzati

Via Bonomea 265, 34136 Trieste, Italy \altaffiliationCurrent address: Theory and Simulation of Materials (THEOS), École Polytechnique Fédérale de Lausanne, 1015 Lausanne, Switzerland

SISSA]SISSA – Scuola Internazionale Superiore di Studi Avanzati

Via Bonomea 265, 34136 Trieste, Italy

CNR-NANO]CNR-NANO, Istituto Nanoscienze, Centro S3

Via Campi 213A, 41125 Modena, Italy

SISSA]SISSA – Scuola Internazionale Superiore di Studi Avanzati

Via Bonomea 265, 34136 Trieste, Italy

1 Introduction

Modeling the optical properties of molecular dyes in solution has long relied on a static picture, whereby thermal fluctuations are thoroughly ignored, either by neglecting the effects of the environment altogether, or by accounting for the dielectric screening of the solvent through some kind of continuum solvation model (CSM),1 such as the polarizable continuum model (PCM),2, 3, 4, 5, 6 the conductor-like screening model (COSMO),7, 8, 9, 10 or the (revised) self-consistent continuum solvation (SCCS) model.11, 12, 13, 14, 15, 16, 17 In both cases, the optical spectra are often obtained at the molecular minimum-energy configuration, as computed in vacuo or in some effective medium. Nevertheless, recent studies have shown the importance of thermal fluctuations on the optical properties of solvated dyes. 18, 19, 20, 21, 56, 22 Thermal fluctuations are accurately described in molecular solvation models (MSMs), in which solvent molecules are treated explicitly, i.e. on the same atomistic ground as the solute (this approach is often referred to as explicit solvent model). Some of the previous attempts to account for thermal fluctuations in the optical properties of solvated molecules were based on molecular dynamics (MD) utilizing classical force fields accurately targeted to the specific system(s) of interest. 18, 21, 22 While this approach is apt to cope with thermal effects on the optical absorption spectra, its practical implementation is hampered by the limited transferability of the available force fields, which would require extensive fitting and validation procedures for each system to be examined. In one of our previous works19 we introduced an MSM approach to simulating the optical properties of solvated dyes where the entire system (solute plus solvent) is treated quantum mechanically and let evolve according to ab initio (AI) molecular dynamics,24, 25 while absorption spectra are computed on the fly using time-dependent density-functional theory (TDDFT),26, 27 and averaged along the trajectory. This naturally accounts for both thermal and dielectric solvation effects on the optical spectra, yet not requiring any unwieldy system-specific fitting or validation. When it comes to simulating large systems comprising hundreds of atoms, however, the deployment of ab initio MSMs is hindered by the large computational cost of hybrid functionals,28, 29, 30, 32 particularly when using plane-wave basis sets, as it is common practice in many applications involving AIMD. In a previous paper56 we proposed to use a GGA exchange-correlation (XC) energy functional33 to generate explicit-solvent AIMD trajectories and a semi-empirical morphing procedure, tailored on the gas-phase B3LYP29 TDDFT spectra of the dye, to correct the GGA spectra evaluated along the the AIMD trajectory. GGA is generally assumed to be accurate enough for the sake of sampling AIMD trajectories, in spite of the minor, but not always negligible, differences existing between the harmonic frequencies computed using different functionals.31 As for the morphing procedure previously proposed, while it proved accurate enough to deal with global properties of the optical spectra, such as the color expressed by a given solution, it may not be so when it comes to predicting individual spectral features.

In this paper we propose an approach to this problem that combines the advantages of CSMs and MSMs, by applying them selectively to different steps of a multimodel simulation protocol, while avoiding or mitigating some of their drawbacks, such as the neglect of thermal fluctuations in CSMs and the high computational cost of MSMs, and without using any semi-empirical procedure. In our approach thermal fluctuations are accounted for by using AIMD in conjunction with MSM for the solvent using GGA energy functional. Optical spectra are then computed on the fly using TDDFT and a hybrid functional such as B3LYP,29 which is well suited for spectroscopic purposes,34 and averaged along the AIMD trajectory. In order to save on computer resources, this step is performed by representing solvent molecules with an effective medium.23

The combination of MSM and CSM methods into the same simulation protocol is the cornerstone of our approach: while CSMs have been extensively tested for individual molecular configurations (e.g. the minimum-energy geometry), the effects of CSMs on configurational averages of the solute spectra have not been addressed before. Hence, our approach is thoroughly tested by comparing its predictions with those of a pure MSM in the paradigmatic case of cyanidin-3-glucoside (cyanin) in its neutral charge state (quinonoidal base, A), which is stable at weakly acidic conditions (see Figure 1). 35, 36 Cyanin features both a $\pi$-aromatic skeleton and lateral hydroxyl terminations, thus being representative of different out-of-plane (hydrophobic) and in-plane (hydrophilic) solute-solvent interactions.37 This comparison, which results in a high level of agreement, provides a solid benchmark of CSMs and demonstrate that, while spectra obtained from CSMs and MSMs may differ considerably for individual configurations along the AIMD trajectory, these differences tend to vanish upon time averaging. The power of our methodology is further and significantly confirmed by comparing its predictions with experiment both at strongly acidic conditions, where the positive flavylium cation (AH+) is stable, and at weakly acidic conditions, where the flavylium cation and the quinonoidal base (A) coexist (see Figure 1). 35, 36

The rest of the paper is organized as follows. In Sec. 2 we introduce our simulation protocol; in Sec. 3 this protocol is validated by a detailed comparison of its predictions for the quinonoidal base A with those of a fully explicit MSM; in Sec. 4 the predictive power of our multimodel approach is demonstrated against experimental absorption spectra, measured at strongly and weakly acidic conditions; Sec. 5 finally contains our conclusions.

2 Methodology

Our multimodel approach consists of the following two-step protocol:

MSM-AIMD – We first perform a tens-of-picoseconds-long AIMD simulation of the dye in solution at ambient conditions, where the solvent is treated explicitly at the MSM level (see Figure 2, left panel). In view of the minor effects of hybrid functionals on the magnitude of thermal fluctations, and of their considerable computational cost, this step can be often conveniently performed at the GGA level of theory. A number of statistically independent molecular configurations (frames) is selected along the trajectory, typically a fraction of ps apart one from the other.

CSM-TDDFT – For each frame thus generated, we then compute the TDDFT spectrum of the system obtained by stripping off the water molecules that surround the solute, and mimicking their effects with a continuum solvation model, such as the SCCS model (see Figure 2, right panel). Due to the poor performance of GGA in computing excited state properties, this step is best performed using a hybrid XC functional such as B3LYP.29, 34 The spectra thus computed for each individual frame are finally time-averaged over the entire trajectory.

2.1 Implementation

The multimodel simulation protocol introduced above has been implemented using the Quantum ESPRESSO suite of computer codes,38 which provides specific modules to perform standard density functional theory (DFT) calculations (pw.x), AIMD simulations (cp.x), and to compute optical spectra within time-dependent density functional perturbation theory (TDDFpT: turboTDDFT).39, 40, 41, 42, 43, 44 The SCCS model used to simulate the effects of dielectric screening on the structure of solvated molecules is implemented in the Environ add-on to Quantum ESPRESSO,15 recently extended to TDDFpT.23

All calculations were performed using the plane-wave (PW) pseudopotential (PP) method, with a tetragonal supercell of dimensions $19.2\times 19.2\times 16.0$ Å3. These dimensions were chosen to accommodate 171 H2O molecules, mimicking the first solvation shell at ambient conditions (Figure 2, left panel). The resulting system consisted of 565 atoms and 1538 valence electrons for neutral cyanin. In the case of flavylium cation the number of electrons is the same, while the structure counts one extra proton (Figure 1). Electrostatic divergences have been removed in this case by introducing a neutralizing uniform background. This simulation setup is adequate to model the effect of the solvent, while keeping the interaction of cyanin with its fictitious periodic images small enough as to be negligible.

We have used the Perdew-Burke-Ernzerhof (PBE)33 GGA XC functional, when doing AIMD simulations, while both PBE and the B3LYP29 functionals have been used for computing optical spectra, for benchmarking and production purposes, respectively. Ultra-soft and norm-conserving PPs from the Quantum ESPRESSO public repository45 were were used in conjunction with the PBE and B3LYP functionals, respectively. Kohn-Sham orbitals have been expanded in PWs up to a kinetic-energy cutoff of 25 and 50 Ry when using ultra-soft and norm-conserving PPs, respectively, and in both cases the electronic charge-density and potentials have been expanded in PWs up to 200 Ry. A reduced cutoff of 50 Ry was used to implement the Fock-exchange operator in the B3LYP functional. 44

In the SCCS model the homogeneous dielectric medium is characterized by the experimental dielectric constant of the solvent, which in the case of water at $T=298$ K equals 78.5 and 1.8 in the static and optical regimes, respectively.46, 47 The solute is hosted by a cavity defined by the molecular ground-state electron charge-density distribution, $n(\mathbf{r})$, such that the dielectric constant goes smoothly from 1 inside (when $n>5\times 10^{-3}\leavevmode\nobreak\ \mathrm{a.u.}$) to the value appropriate to the solvent outside (when $n<10^{-4}\leavevmode\nobreak\ \mathrm{a.u.}$). 15, 23 The resulting scheme is illustrated in Figure 2 (right panel). Only electrostatic screening is considered in our protocol, while all non-electrostatic effects (such as Pauli repulsion, dispersion, cavitation, etc.) are consistently disregarded.

AIMD was performed using the Car-Parrinello extended Lagrangian formalism,24, 48 with a fictitious electronic mass $\mu=\mathrm{600\leavevmode\nobreak\ a.u.}$ and a time step $\Delta t=\mathrm{0.12\leavevmode\nobreak\ fs}$. In order to facilitate the approach to equilibrium and enhance dynamical sampling, all the atomic masses were set to a same value of $\mathrm{M=12\leavevmode\nobreak\ a.m.u.}$ Dispersion forces have been accounted for by using the semi-empirical Grimme’s DFT-D2 van der-Waals interaction correction.49, 50 Our AIMD run consisted of an equilibration phase of $\sim 2.5\leavevmode\nobreak\ \mathrm{ps}$ at $300$ K using a Nosé-Hoover thermostat,51, 52, 53 followed by a 50 ps production trajectory. The AIMD trajectory thus generated was sampled every $\mathrm{0.5\leavevmode\nobreak\ ps}$,111We estimate that the auto-correlation time of the HOMO-LUMO gap is smaller than 100 fs, while the typical auto-correlation times of solvent-solute hydrogen bonds is of the order of the ps. resulting in 100 statistically independent frames for which the absorption spectra have been computed by either keeping water molecules at their instantaneous positions (MSM), or by representing them with an effective continuous medium (CSM). The calculation of the absorption spectrum at every frame has been performed in two steps: i) we have first computed the ground-state electronic structure, using the pw.x code; ii) then, we have computed excitation energies and oscillator strengths by solving Casida’s equation54 with a Davidson-like algorithm55, 44 implemented in the latest release of the turboTDDFT code.43, 44 Absorption spectra have been convoluted with a Lorentzian broadening function with a full width at half maximum of $\approx 0.1$ eV. The RGB representation of the color expressed by solvated dyes has been obtained by using the tristimulus colorimetry theory, as described e.g. in Refs. 19 and 56.

For the sake of comparison, we have also calculated the absorption spectra of the molecule in vacuo, at the various geometries generated by the AIMD trajectory upon removal of the water molecules (we refer to the resulting spectra as the dehydrated spectra). We stress that no structure relaxation has been performed for these structures, which were thus kept fixed at the instantaneous geometries generated by AIMD, so as to account for the effects of thermal vibrations.

3 Validating the multimodel approach

In this section the quality of our multimodel approach is assessed against the predictions of a full MSM simulation of the absorption spectra. For the sake of validation, we consider the case of the quinonoidal base A. The comparison is done at the same level of theory. The size of the explicit model used here (171 water molecules, see above) would make B3LYP calculations extremely expensive. On the other hand, the computation of absorption spectra for this model can be afforded using the PBE functional, which, although inadequate for quantitative predictions, can indeed capture the qualitative features of the optical spectra of anthocyanins.57 In fact, while GGA XC functionals underestimate the HOMO-LUMO gap and other low-lying transition energies, the relative order of the the first few of them, as well as the character of the molecular orbitals participating in them are correctly predicted. For this reason our assessment was performed using the PBE functional, while a B3LYP functional will be used for the comparison with experiments (see Sec. 4), for which no full MSM simulation was attempted.

We start by comparing the lowest energy (i.e. longest-wavelength) transitions of the molecule in the gas phase, as computed at the minimum-energy geometry configuration in vacuo, with that resulting from CSM, as computed at the SCCS minimum-energy configuration. The two spectra (Figure 3, upper panel) are characterized by an intense peak at $\lambda_{\mathrm{max}}\approx 570\leavevmode\nobreak\ \mathrm{nm}$.23 Accounting for dielectric screening results in a blue shift of the first absorption peak by a few nanometers only and to an increase of its intensity.

In order to evaluate the impact of thermal fluctuations independently of solvent polarization effects, in Figure 3 (lower panel) we report the dehydrated spectrum (i.e. obtained by averaging the spectra computed in vacuo for the molecular geometries sampled by the AIMD trajectory), featuring a large broadening and a red shift of the first absorption peak by $\Delta\lambda_{\mathrm{max}}\approx 50\leavevmode\nobreak\ \mathrm{nm}$. An explicit account of both dielectric-screening and thermal effects is achieved in both the MSM and multimodel approaches, as illustrated in Figure  3, resulting in a markedly different behavior from that obtained by neglecting either one of these effects. The time-averaged spectra from MSM and multimodel simulations agree remarkably well, with minor residual differences consisting in a deviation in $\lambda_{\mathrm{max}}$ of a few nm, and in a more pronounced broadening predicted by MSM. The latter effect is likely due to the larger fluctuations of the electrostatic potential induced by the molecular nature of the solvent and/or to the neglect of non-electrostatic solute-solvent interactions in our implementation of the SCCS model.

Figure 4 reports the MSM, CSM, and dehydrated spectra computed at a few individual and non-consecutive frames, representative of the full AIMD trajectory. Comparing the spectra obtained at different levels of approximation for different frames does not reveal a consistent picture. Sometimes different approximations result in very similar spectra (Frame A); in the majority of cases, however, the value of $\lambda_{\mathrm{max}}$ predicted by different solvation models do not coincide (Frames B – E) and, moreover, the shape of the spectrum can differ considerably, featuring one or two intense peaks in the range 500 – 800 nm (Frames B and D). The difference between MSM and CSM spectra for a same configuration of the solute is obviously a manifestation of the discrete (molecular) nature of the solvent. In the present case, inspection of the microscopic configurations of water molecules around the solute reveals a noticeable correlation between the number of hydrogen bonds formed by water molecules and the most acidic O atom in position 4’ (see Fig. 1) and the difference between the MSM and CSM wavelengths of the fundamental absorption peak, $\lambda_{\mathrm{max}}$.

It is interesting to compare the absolute intensities of the spectra obtained from different solvation models: the resulting trend $I_{\mathrm{MSM}}>I_{\mathrm{CSM}}>I_{\mathrm{deh}}$ can be ascribed to the larger delocalization of molecular orbitals in vacuo than in solution, which implies a smaller excitation oscillator strength. Molecular orbitals are more delocalized in CSM than in MSM, because the latter model explicitly accounts for Pauli repulsion from the solvent molecules, whereas the former does not (at least in our implementation).

Checking the convergence of our time-averaged absorption spectra requires some care. In Figure 5 we display the absorption spectra of cyanin computed in MSM and CSM using the PBE functional, and averaged over an increasing number of consecutive frames, evenly spaced by 0.5 ps. We see that up to $\approx 60$ independent frames are needed to achieve a reasonable convergence. When not otherwise specified, the spectra presented in the text are obtained by averaging over 100 configurations. The convergence of the time-averaged spectrum in CSM using B3LYP XC functional follows the same trend.

4 Predicting experimental spectra

Figure 6 (upper panel) shows the spectra predicted by our multimodel approach for the solvated cyanin, in both the flavylium cation (pink), and quinonoidal neutral base (purple) charge states, using the B3LYP XC functional. At very low pH (<3), the flavylium cation, is the only stable species present in solution, and our theoretical predictions for its spectrum can thus be directly compared with experiment (Figure 6, middle panel). The agreement with the experiments is very good as concerns both the overall shape of the spectrum, including the position of the maximum-absorption peak ($\lambda_{\mathrm{max}}^{\mathrm{theor}}=520$ nm, $\Delta\lambda_{\mathrm{max}}\approx 12$ nm), and the hue of its simulated color. Although our simulations predict a secondary absorption peak at shorter wave-length ($\lambda\approx 425\leavevmode\nobreak\ \mathrm{nm}$) that is not observed experimentally, a shoulder is clearly visible in the experimental spectrum at this wave-length.58

As the acidity of the solution decreases, different molecular species coexist: i) deprotonation of the phenyl ring can take place giving rise to the quinonoidal base; ii) hydration of the quinonoidal base generates several optically inactive, colorless species (e.g. chalcones).35 The pH-dependent molar fractions of these species have been experimentally determined in Ref. 59. At $\mathrm{pH}=4$ the relative concentrations of the various species involved in the chemical equilibrium were found out to be: $\mathrm{\left[AH^{+}\right]=3\%}$, $\mathrm{\left[A\right]=5\%}$, and $\mathrm{\left[CL\right]=92\%}$, where [CL] indicates the molar fraction of the ensemble of the colorless components. The high percentage of colorless components bleaches the color of the overall solution. A residual optical activity, however, still persists, and our predictions for the optical properties of the quinonoidal base can be validated by comparing the weighted sum of the spectra computed for the flavylium cation and quinonoidal base with the measured absorptivity. Our results are shown in Figure 6 (lower panel), which displays once again a very good agreement with experiments:58 $\lambda_{\mathrm{max}}^{\mathrm{theor}}=529$ nm, $\Delta\lambda_{\mathrm{max}}\approx 8$ nm (actually even better than at very low pH, probably as a result of a fortuitous cancellation of errors).

5 Conclusions

In this paper we have introduced a new approach to simulate the optical properties of complex dyes in solution in the visible range, based on a multimodel protocol to account for the effects of the solute-solvent interactions (such as thermal vibrations and dielectric screening) on the absorption spectra of molecular dyes. Our methodology builds on previous attempts to simulate thermal fluctuation effects using classical18, 21, 22 or ab initio19 molecular dynamics and TDDFT to compute the optical spectra for individual molecular configurations thus generated, and improves upon them by avoiding to use classical force fields, one one hand, and by optimally combining explicit and implicit solvation models on the other hand. The validation of the newly introduced protocol, performed on cyanidin-3-glucoside in its neutral (quinonoidal base) charge state at room temperature, sheds light onto the validity of implicit (continuum) solvation models of the structural and spectroscopic properties of molecular dyes in solution, and results in a remarkable agreement with the experiments in the strong-to-weak acidity range.

Our developments open the way to the simulation of the optical properties of complex solvated dyes in more realistic conditions, and with a higher accuracy than it has been possible so far. Our approach, which is based on extensive sampling of the optical properties of the solute corresponding to different molecular configurations occurring at equilibrium, is being further enhanced by accounting for low-frequency conformational fluctuations occurring over widely varying time scales, 60 by developing new “learn-on-the-fly” techniques to save on the time-consuming TDDFT calculations over molecular dynamics trajectories,61 and by combining it with a QM/MM approach to solute-solvent interactions.

{acknowledgement}

We are grateful to Rebecca Robbins and Tom Collins for illuminating discussions and for providing us with their unpublished data of Ref. 58. We also thank Xiaochuan Ge and Oliviero Andreussi for cooperating with us during some of the preliminary phases of the present work, and Carlo Cavazzoni and Fabio Affinito for technical support at CINECA. This work was partially supported by Mars Chocolate North America LLC and by the European Union through the MaX Centre of Excellence (Grant No. 676598). Computer resources were partially provided by the CINECA Supercomputing Center (Italy) at their Fermi BG/Q Tier-0 and Galileo Tier-1 machines through grants PRACE No. 2013081532 (Chromatology), IscrC_MULAN, IscrC_SCCS-H, and IscrC_SCCS-H2.