Elizete Ventura*a,
Silmar Andrade do Monte
*a,
Mariana T. do Casal
b,
Max Pinheiro Jr
b,
Josene Maria Toldo
b and
Mario Barbatti
*bc
aUniversidade Federal da Paraíba, 58059-900, João Pessoa-PB, Brazil. E-mail: elizete@quimica.ufpb.br; silmar@quimica.ufpb.br
bAix Marseille University, CNRS, ICR, Marseille, France. E-mail: mario.barbatti@univ-amu.fr; Web: http://barbatti.org
cInstitut Universitaire de France, 75231 Paris, France
First published on 5th April 2022
The heating of a chromophore due to internal conversion and its cooling down due to energy dissipation to the solvent are crucial phenomena to characterize molecular photoprocesses. In this work, we simulated the ab initio nonadiabatic dynamics of cytosine, a prototypical chromophore undergoing ultrafast internal conversion, in three solvents—argon matrix, benzene, and water—spanning an extensive range of interactions. We implemented an analytical energy-transfer model to analyze these data and extract heating and cooling times. The model accounts for nonadiabatic effects, and excited- and ground-state energy transfer, and can analyze data from any dataset containing kinetic energy as a function of time. Cytosine heats up in the subpicosecond scale and cools down within 25, 4, and 1.3 ps in argon, benzene, and water, respectively. The time constants reveal that a significant fraction of the benzene and water heating occurs while cytosine is still electronically excited.
Our present work aims to quantify the time scales for both processes, including excited-state energy transfer and nonadiabatic events. Taking cytosine as a prototypic chromophore undergoing ultrafast internal conversion, we simulated its nonadiabatic dynamics in an argon matrix, benzene, and water (Fig. 1); thus, we gauged the energy transfer from weak van der Waals’ interactions, through medium intensity π–π interactions, to strong hydrogen bonds. Then, we proposed an analytical model to quantify the transfer rates, providing a protocol that can be used with data from any other dynamics simulations.
![]() | ||
Fig. 1 Systems studied in this work: clusters of (a) cytosine and Ar atoms; (b) cytosine and benzene molecules, and (c) cytosine and water molecules. |
Not unexpectedly, we are not the first group to address the fundamental physical–chemical problem of heat transfer between a molecule and its solvent.1–10 Nevertheless, studies that describe solute heating and cooling in solution after an initial electronic excitation are rare.1 Most theoretical models rely on classical dynamics and only describe solute cooling. The cooling process is investigated in these models, supposing that the solute initially has high kinetic energy2–6 or high potential energy.7,8 In some models used to study intramolecular vibrational relaxation, the photon energy is assumed to be instantaneously converted into the vibrational energy of a given moiety.9,10 Recently, Balevičius Jr et al.1 proposed an approach in which a simplified Hamiltonian is used to describe both heating and cooling processes. Their model, however, does not take into account nonadiabatic and leaking effects. This latter includes the possibility that the initial excess (electronic) energy relaxes simultaneously to the internal degrees of freedom of the solute and solvent.
Our motivation to develop this project has been connected to our latest works on the development of molecular heaters for agriculture in extreme weather.11 Nevertheless, the range of interest of our results extends much beyond this topic. Recent technological applications profit from chromophore-solvent heat transfer, such as photothermal therapy (where the heat generated by the photoexcited chromophore is used to kill cancer cells and pathogens)12,13 and molecular solar thermal energy storage (where a high-energy metastable isomer can release the stored energy to the environment as heat).14,15 In contrast, heat release to the environment is undesirable in many other applications, such as singlet oxygen generation in photodynamic therapy,14,16,17 data storage,18 and solar light-harvesting applications.19,20 Additionally, the study of chromophore-solvent heat transfer is essential to understand the reactivity of molecules in solution, as the rates of these reactions can be strongly dependent upon the internal energy of the reactants and hence upon the energy flow between solute and solvent molecules.21 Indeed, in the presence of competing processes, the chromophore-solvent heat-transfer rates can be decisive, controlling the reaction yields.
DNA nucleobases and nucleosides are excellent models for studying the vibrational energy flow from photoexcited molecules to their surroundings because their ultrafast and complete internal conversion to S0.22–24 We have chosen cytosine as a prototypical chromophore because its internal conversion is well-characterized theoretically25–40 and experimentally.31,33,39–44 After light absorption, its excited-state dynamics is governed by a branching between competing decay channels, associated with three conical intersections, ππ*/S0 (ethylenic), nOπ*/S0, and nNπ*/S0. Although there is some debate about which one of the conical intersections dominates the nonradiative decay in the gas phase, several theoretical25,37,45,46 and experimental works43,47,48 agree on reporting ultrafast time constants (τ1 = 0.01–0.2 ps and τ2 = 0.5–2 ps), primarily associated with the decay through one of the nπ*/S0 conical intersections. The ultrafast nonradiative decay is connected with its photostability and effectively dissipates the harmful electronic excitation into heat.34
Although there are systems that show intermolecular vibrational relaxation mainly in the excited states (like perylene in ketone solvents49), the ultrafast internal conversion of cytosine implies a competition between relaxation from the excited and ground states. From the theoretical point of view, describing such a competition requires multiconfigurational methods to account for the nonadiabatic dynamics to the ground state.50 The approach used in this work (employing surface hopping simulations coupled with hybrid CASSCF/molecular-mechanics electronic structure51,52) can describe the relaxation independently of its time scale and source state (ground or excited). Thus, it can deal with the nonadiabatic relaxation of the initial electronic energy to the internal degrees of freedom of the solute and solvent of a realistic system, features that are not considered in most of previous approaches.
![]() | (1) |
![]() | (2) |
![]() | (3) |
The different pathways described by each pi connect the first excited state to the ground state through different time constants τci. In eqn (1), all three solvents could be described by two reaction pathways (Np = 2).
The term in the parentheses on the right side of eqn (1) transfers energy to the solvent (for Ec > NcEs/Ns). It is the difference between the chromophore's kinetic energy and the solvent's mean kinetic energy over the same number of atoms. (This term is analogous to what we expect from thermal contact between a hot and a cold body.) τC is the chromophore's cooling time constant, while Nc and Ns are the numbers of atoms in the chromophore and solvent, respectively.
In principle, Es should evolve according to a differential equation coupled to eqn (1). However, if we assume that there is little energy back transfer from the solvent to the chromophore, we can first find out the functional dependence of Es on time and use it to solve eqn (1). Our simulations show that
![]() | (4) |
![]() | (5) |
Employing eqn (3) and (4), the solution of eqn (1) is
![]() | (6) |
![]() | (7) |
![]() | (8) |
The coefficients Ai can either be positive or negative, reflecting the balance of heat gain or loss at that time constant. Assuming thermal equilibrium between chromophore and solvent at t → ∞ and that there is a harmonic virial partition between kinetic and potential energies, the asymptotic kinetic energy values of the chromophore and solvent are connected to the initial values through
![]() | (9) |
![]() | (10) |
In these equations, hν is the photon energy, E0c = Ec(0), and E0s = Es(0). A general discussion about these equilibrium values is in the ESI† Fig. S1.
Returning to eqn (6), C depends on Ai and Bn due to the boundary values of Ec at t = 0 and t → ∞:
![]() | (11) |
To use the model, we proceed through the following steps:
1. We fit the ground-state population from nonadiabatic dynamics simulations with eqn (2) to obtain ai and τci;
2. We use the kinetic energy of the solvent, also coming from dynamics simulations, to determine E∞s with eqn (10) and get κsn and τsn through a fitting procedure;
3. We use the kinetic energy of the chromophore to get E∞c (eqn (9)) and fit these data with eqn (6) to determine κci and τC. These parameters are constrained to be ≥0.
This procedure allows determining the time to heat the chromophore due to internal conversion and the time to cool it down due to energy transfer to the solvent. The heating time is given by the mean excited-state lifetime
![]() | (12) |
Each τsn in eqn (4) is a time constant measuring a single channel of solvent heating. In turn, τC is a single time constant describing the overall solute cooling. In a fully coupled model, these constants balance, and a common time constant given by the harmonic mean of all these time constants emerges. Our model, however, neglects this back coupling by proposing a solvent kinetic energy dependence (eqn (4)) independent of the solute. Thus, the complementarity between time constants is lost. If the solvent heats through a single channel and there is no heat via internal conversion, we could impose this complementarity as an additional constraint. However, we do not see how to do it when multiple channels are present. We will see in the results that neglecting the back couplings is justifiable in all three cases we examined. Thus, we focused on τC time constant as the measure of the solute cooling time.
The solvent was treated via molecular mechanics (MM). Three solvation schemes were considered: (1) cytosine in a 20 Å sphere containing 684 Ar atoms at 10 K; (2) cytosine in a 19.2 Å sphere with 200 benzene molecules at 298 K; (3) cytosine in a 12.9 Å sphere with 300 water molecules at 298 K. The OPLS/AA force field was used for Ar, cytosine, and benzene. For Ar, standard OPLS/AA55 parameters were used, while for cytosine and benzene, they were obtained using the LigParGen web server.56–58 The TIP3P59 force field was used for water. The solute–solvent interaction was computed through QM/MM in an electrostatic embedding.60 Cytosine was in the QM region and the solvent in the MM region.
Solute–solvent (microcanonical) nonadiabatic dynamics simulations were done with QM/MM surface hopping. The hopping probabilities were computed with the decoherence-corrected61 fewest-switches surface hopping62 (DC-FSSH). The quantum integration is done with 0.025 fs using interpolated electronic quantities between classical steps of 0.5 fs. The initial conditions were generated with the hybrid Wigner-thermal protocol proposed in ref. 52. For the simulations in argon and benzene, we ran 50 trajectories each. In water, we ran 75 trajectories. The initial state was distributed over S1 to S3 according to the transition probabilities in the 5.25 ± 0.25 eV energy window.
All CASSCF calculations were done with the COLUMBUS program.63–66 The MM calculations were done using TINKER software.67 Surface hopping dynamics was performed using the NEWTON-X software68 interfaced with COLUMBUS and TINKER.
Table 1 shows the weights (ai), time constants τci (see eqn (3)), and cytosine's average excited-state lifetime in the three solvents and in the gas phase (from ref. 25). These results are obtained by fitting the fraction of trajectories in the ground state as a function of time with eqn (2) up to 1.0 ps in the gas phase, 3 ps in argon, 1.5 ps in benzene, and 1.5 ps in water (see Fig. 2).
![]() | ||
Fig. 2 Ground-state population as a function of time simulated with surface hopping for cytosine in argon (top), benzene (middle), and water (bottom). The dashed lines are the fitting function from eqn (2). |
The QM/MM mean kinetic energies of cytosine Ec(t) and each solvent Es(t) are shown in Fig. 3. The results of the energy-transfer model (eqn (4) and (6)) are plotted too (Nc = 13 and hν = 5.25 eV). The agreement between them is outstanding. The bottom graphs in Fig. 3 show the long-timescale behavior of the energy-transfer model. The energy distribution between translational and internal degrees is examined in ESI† Fig. S4.
![]() | ||
Fig. 3 Time dependence of the mean kinetic energies of cytosine and solvents. The upper graphs show the results from the dynamics averaged over the trajectories. The dashed and dotted lines are the energy-transfer functions of eqn (6) for cytosine and eqn (4) for the solvent. These same functions are shown in the bottom graphs but over an extended time. The kinetic energy of the solvent is scaled by Nc/Ns, like in eqn (1). |
The parameters of the energy-transfer model describing the solvent's kinetic energy Es(t) are given in Table 2. Argon's kinetic energy mainly grows with a single time constant of 15.5 ps. It is much longer than the 0.64 ps excited-state lifetime (Table 1). This result implies that the energy transfer to the solvent happens mostly after cytosine returns to the ground state. In the case of water, the situation is the opposite: water's kinetic energy grows within 0.15 ps, which is much shorter than the excited-state lifetimes (0.62 ps), meaning that the energy transfer starts while cytosine is still excited. Benzene kinetic energy grows in the first 0.5 ps above the equilibrium value and converges to it from above (Fig. 3, middle-bottom). We have tracked this behavior to a mismatch between the MM-equilibrated initial conditions and the QM/MM-level dynamics. Although this artifact rendered too short τs1 (0.06 ps in Table 2), it does not impact our main results, the description of Ec(t), because Ec(t) ≫ EsNc/Ns in eqn (1).
Solvent | |||
---|---|---|---|
Argon | Benzene | Water | |
N s | 684 | 2400 | 900 |
E 0s (eV) | 0.89 | 92.03 | 34.66 |
E ∞s (eV) | 3.46 | 94.64 | 37.25 |
κ s1 (eV) | −2.576 | −2.61 | −2.59 |
τ s1 (ps) | 15.49 | 0.06 | 0.15 |
The energy-transfer parameters for cytosine are given in Table 3. In the three solvents, B1 ≪ Ai, which is consequence of the much larger kinetic energy (per atom) in cytosine than in the solvent. It also implies a low level of back energy transfer, as we assumed in the model. In water, we additionally observe A2 = 0, meaning that the cytosine's heating through the second decay pathway is canceled out by its simultaneous energy transfer to the solvent.
Solvent | |||
---|---|---|---|
Argon | Benzene | Water | |
E 0c (eV) | 1.46 | 1.46 | 1.43 |
E ∞c (eV) | 1.50 | 1.47 | 1.47 |
A 1 (eV) | −0.97 | −1.37 | −1.72 |
A 2 (eV) | −1.55 | −1.74 | 0.00 |
B 1 (eV) | 0.08 | 0.00 | 0.005 |
C (eV) | 2.40 | 310 | 1.68 |
τ C (ps) | 24.9 | 3.7 | 1.3 |
The cytosine cooling time, τC, is also given in Table 3. As expected, the energy transfer in argon is the slowest, taking about 25 ps. This time drops to about 4 ps in benzene and is merely 1.3 ps in water.
Our estimate for the cooling time of cytosine in benzene is 4 ps (Table 3). For comparison, after excitation at 261 nm, the cooling of coumarin in cyclohexane is about 10 ps, as measured with time-resolved fluorescence spectroscopy.69 The chromophore-solvent interaction between cytosine and benzene should be stronger than that between coumarin and cyclohexane (as the π–π interaction is absent in the latter), explaining the slightly shorter time for cytosine.
Zhang et al.,70 using UV-pump/broadband-mid-IR-probe transient absorption spectra, associated the cooling time of hot purine derivatives to the directionality of their H-bonds, showing that molecules that have more N–H bonds have shorter cooling times. They proposed that N–H(D) bond donation is responsible for rapid energy disposal to water via direct coupling of high-frequency solute–solvent modes. Thus, the cooling times decrease in the order caffeine (7.7 ps) > theophylline (5.1 ps) > hypoxanthine (3 ps) due to the larger number of N–H bonds.70 Our value of 1.3 ps obtained for cytosine in water (Table 3) fits very well in this sequence, as our chromophore makes three N–H bonds, while hypoxanthine has two. Our cooling time is consistent with the vibrational cooling of 2 to 3 ps of photoacid pyranine71 and 2.4 ps of 9-methyl-adenine,72 both in water. Foremost, it is consistent with the vibrational cooling of cytosine in water, which is reported as ∼2.9 ps23,31 and ∼4.0 ps.73 Indeed, it has been shown that hydrogen bonds between nucleobase monomers and solvent strongly enhance vibrational cooling.22,31 The importance of H-bonds for a faster solute–solvent vibrational energy relaxation has also been pointed out by Pigliucci et al.,74 who identified that such process is enhanced for solutes (substituted perylenes) bonded to alcohols.
Cytosine's kinetic energy reaches a maximum after 2.38 ps in argon, 0.95 ps in benzene, and 0.16 ps in water, as we can see in Fig. 3. This time represents a transient equilibrium between the energy that the internal conversion transfers to cytosine's vibrational modes and the energy cytosine transfers to the solvent. Interestingly, this transient equilibrium in water happens before the internal conversion (0.62 ps), meaning cytosine starts to cool down while still photoexcited. The outcome for water is consistent with the findings of ref. 74, which concluded that both intra- and intermolecular (solute–solvent) vibrational energy redistribution processes occur, at least partially, on similar timescales.
We do not expect that the modest computational levels employed in this work will be the last word in the description of cytosine-solvent transfer times, although the analysis above shows a semi-quantitative agreement with previous results. Nevertheless, a remarkable feature of our model is that it is independent of the Hamiltonian, and it can be employed with more advanced methods,75,76 as long as nonadiabatic dynamics is affordable.
We have applied this energy-transfer model to a chromophore embedded into three types of solvents, argon, benzene, and water, spanning an extensive range of interactions. Due to its complete and ultrafast (sub-picosecond) internal conversion, we adopted cytosine as a prototypical chromophore. Using data from surface hopping dynamics with QM/MM, our model predicts that cytosine cools down within 25 ps in argon, 4 ps in benzene, and 1.3 ps water.
Our initial goal in this work was to estimate the chromophore-to-solvent heat transfer times. These results for cytosine in argon, benzene, and water may be taken as a qualitative indication of the orders of magnitude of the cooling time of other organic chromophores in similar solvents. Nevertheless, the development of an energy-transfer model that can be directly employed with results from any nonadiabatic dynamics simulations goes much beyond that initial goal and delivers a new protocol to analyze energy transfer in complex environments. Currently the model is being extended to also take into account the heat transfer between strongly coupled subsystems (e.g. through a covalent bond), as in the case of nucleosides.
Footnote |
† Electronic supplementary information (ESI) available: Definition of equilibrium energy, computational details, description of the nonadiabatic dynamics for cytosine, decomposition of internal and translational energies. See DOI: 10.1039/d2cp00686c |
This journal is © the Owner Societies 2022 |