Modeling the heating and cooling of a chromophore after photoexcitation

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.


S1. Kinetic energy partition
When t →, the photon energy h should distribute among all degrees of freedom. If we assume equipartition, each degree should receive ( ) The photon contribution for the equilibrium kinetic energies of the chromophore is ( ) ( ) The equilibrium kinetic energy of the chromophore ( ) xc = and solvent ( ) In all three cases treated in our paper-argon matrix, liquid benzene, and water-2 Thus,

A. QM Calculations
The ground state geometry of cytosine was optimized at Complete Active Space Self-consistent Field (CASSCF) level. The complete active space consists of fourteen electrons in ten orbitals, i.e., a CAS (14,10). These ten orbitals at the ground state minimum geometry are three n, four π, and three π* orbitals, described in detail in ref. 1 The first four singlet states were averaged with equal weights. A normal-mode analysis was also performed to confirm that the optimized geometry is a minimum. The geometry optimization and normal-mode calculation were performed using analytical gradient techniques. [2][3][4][5] Frequency calculations were performed via finite differences from analytic gradients. All CASSCF calculations were done with the COLUMBUS program. [6][7][8][9] The atomic orbitals (AO) integrals and their gradient integrals used by COLUMBUS were computed with the corresponding modules taken from the DALTON program. 10 Due to the extensive computational resources required by the dynamic calculations (up to 3ps) repeated over three different solvents, we adopted the 3-21G basis set. 11 Dynamics simulations with this basis set are in semi-quantitative agreement with the results using the 6-31G* basis set reported in ref. 1 and present a significant computational cost reduction. Moreover, since our goal is not to discuss the internal conversion of cytosine but the energy transfer to the solvent, we only need a computational level that can successfully describe the ultrafast dynamics to S0.

B. MM and QM/MM calculations
To build the cluster with the cytosine molecule inside an Ar matrix, a similar procedure as described in ref. 12 was used. First, the Ar matrix was built by increasing the unit cell of the experimental X-ray structure 13 to dimensions 50 x 50 x 50 Å (yielding 1372 Ar atoms), using the vegaZZ software. 14 The same program was used to insert the cytosine molecule inside the cavity. To ensure that the cavity was large enough to accommodate the solute, the default maximum overlap (0.80Å) between solute and neighboring Ar atoms in the matrix was applied, leading to a removal of four Ar atoms from the center of the crystal matrix. Then, a sphere of 20 Å containing the cytosine and 684 Ar atoms was built for the subsequent equilibration (described later). This procedure was carried out only once for one of the selected Wigner structures of cytosine. The remaining clusters were generated by replacing the cytosine structure with each of the other selected cytosine geometries generated by the Wigner distribution.
Spherical clusters of cytosine in benzene and water were built using the PACKMOL package, 15 with radii of 19.2 Å and 12.9 Å, respectively. The number of solvent molecules (200 and 300 for benzene and water, respectively) was chosen based on their densities at 298 K, 0.876 (benzene) and 0.997 (water) g/cm 3 . 16 As for the Ar cluster, such procedure was done only once, using one of the Wigner structures generated for cytosine (see Section C). For the remaining ones, the cytosine structure was replaced accordingly.
All MM calculations were done using TINKER software. 17 Once the cluster structures were obtained, they were thermalized at the MM level using the NVT ensemble, keeping the cytosine frozen. A temperature of 10 K was used for the Ar cluster, while 298 K was used for benzene and water clusters. The equilibration times varied from 5 to 150 ps, depending on the solvent. The Ar clusters required larger equilibration times. The "wall" keyword in TINKER was used during the equilibration to maintain droplet boundary conditions, preventing solvent atoms from moving outside the defined spherical radius. The OPLS/AA force field was used for Ar, cytosine, and benzene. For Ar, standard OPLS/AA 18 parameters were used, while for cytosine and benzene, they were obtained using the LigParGen web server. [19][20][21] The TIP3P 22 force field was used for water. Once the clusters were equilibrated, the initial coordinates of solvent molecules and solute (this latter generated by the Wigner distribution) needed for nonadiabatic QM/MM dynamics were defined. The initial velocities, also required for the initial conditions, were taken from the Wigner distribution for cytosine and MM equilibration for solvent molecules. This procedure, proposed in ref. 23 , prevents the nuclei of the solute from having energies below the zero point but ensures that the solvent is at the correct temperature.
The solute-solvent interaction was computed through QM/MM in an electrostatic embedding. 24 Cytosine was in the QM and the solvent in the MM regions.

C. Initial conditions for dynamics
Initial geometries and velocities for the QM region (cytosine) were generated by a

D. Nonadiabatic dynamics
Nonadiabatic dynamics was simulated with QM/MM surface hopping. The hopping probabilities were computed with the decoherence-corrected 25 fewest-switches surface hopping 26 (DC-FSSH) algorithm. The quantum integration was done with 0.025 fs using interpolated electronic quantities between classical steps of 0.5 fs. Time-dependent electronic coefficients were adjusted for decoherence with the simplified decay of mixing method. 25 The decoherence parameter was set to 0.1 au. In the case of hopping, the momentum excess is adjusted in the direction of the nonadiabatic coupling vectors, while in the case of frustrated hopping, the momentum is kept constant. Additional details concerning the integration of classical and quantum equations are given in ref 27 . Whenever a given trajectory remained for at least 50 fs in the ground state, it was restarted adiabatically in this state. Anelastic-collision spherical boundary was applied to keep the density constant. The sphere is centered at cytosine's center of mass with the cluster radius.
Surface hopping dynamics was performed using the NEWTON-X software 28  As discussed in Section C, we have yielded 79 initial conditions, 2 starting in S1, 47 in S2, and 30 in S3. For the simulations in argon, we ran 50 trajectories (1 starting in S1, 30 in S2, and 19 in S3). In benzene, we also ran 50 trajectories with the same distribution of initial states. In water, we ran 75 trajectories (2 starting in S1, 45 in S2, and 28 in S3). The number of trajectories in each initial state is summarized in Table 1.
In argon, 30 trajectories run at least 1.8 ps, and 18 ended at 3 ps. In benzene and water, 30 trajectories of each solvent ran at least 1.5 ps. In argon, the heat-transfer model was fitted with data up to 3 ps. In benzene and water, the fitting was done with data up to 1.5 ps.

S3. Dynamics of cytosine
As expected, the time constants and the excited-state lifetime of cytosine in argon matrix and benzene are close to that obtained in the gas phase. 1 However, in water, the corresponding time constants are considerably larger, although the average excited-state lifetime (0.62ps) is also close to those obtained in Ar and benzene. According to previous surface hopping simulations in the gas phase, 1 the larger time constant could be interpreted as a decay channel mainly through a semi-planar conical intersection along the nOπ* surface. In fact, from the fraction of trajectories that ended up adiabatically in S0, i.e., the ones which were restarted as adiabatic dynamics in the S0 after the S1/S0 crossing (90, 78, and 76% in Ar, benzene, and water, respectively), 97% (Ar), 85% (benzene), and 100 % (water) decayed through this semi-planar conical intersection.
The weight and the time constant associated with the ultrafast decay channel (< 50 fs) increase in the following order: argon < benzene < water. These two results explain the average excited-state lifetime of cytosine in water being close to those in the other two solvents, despite its significantly larger 2 c  value. The second most important conical intersection associated with the S1→S0 decay is the C6-puckered, 1 also reported as an ethylenic conical intersection, which involves the decay of the * state. 29 Computational simulations including explicit water molecules support the involvement of the n* state. 31,33 Although the ultrafast time constants of cytosine in water are often linked to the decay of * state, the agreement between experimental and our theoretical results is reasonable, considering the electronic-structure level and the number of trajectories in our simulations. Therefore, because our simulations do not include explicit water molecules in the QM part, it is natural that our results resemble the ones in the gas phase, 1 since quantum-mechanical interactions, such as water-chromophore electron transfer, 33 are not captured at this level of theory. Nevertheless, the strength of the cytosine-water hydrogen bonds is correctly predicted by our QM/MM partition, which is enough to ensure the energy flow between the subsystems. This region includes trajectories up to 3.0 ps for cytosine in Ar and 1.5 ps for cytosine in benzene and water.
It is worth mentioning that a fraction of trajectories (6 to 20%, depending on the system) remained in the excited states, similar to earlier surface hopping simulation in the gas phase. 1 This allows an alternative interpretation, using a tri-exponential fitting, where part of the population would decay with longer time constants.

S4. Internal and translational energies
We  The initial reduction of the ratios shown in Figure 2 is an artifact due to the MM equilibration performed on the clusters before they were submitted to the nonadiabatic QM/MM dynamics. The initial decrease can be due to a relaxation of the system under the influence of the new QM/MM environment. It is expected that with an initial equilibration at the QM/MM level, 12 the ratio should show a smaller variation.