We systematically applied excited-state normal mode analysis to investigate and compare the relaxation and internal conversion dynamics of a free-base porphyrin (BP) with those of a novel functional porphyrin (FP) derivative. We discuss the strengths and limitations of this method and employ it to predict very different dynamical behaviors of the two compounds and to clarify the role of high reorganization energy modes in driving the system toward critical regions of the potential energy landscape. We identify the modes of vibrations along which the energy gap between two excited-state potential energy surfaces within the Q band manifold may vanish and find that the excess energy to reach this “touching” region is significantly reduced in the case of FP (0.16 eV) as compared to the one calculated for BP (0.92 eV). Our findings establish a link between the chemical functionalization and the electronic and vibrational structure that can be exploited to control the internal conversion pathways in a systematic way.

## I. INTRODUCTION

Studying the synthesis, photo-physics, and photo-chemistry of porphyrin derivatives is a long-standing research topic, which has gained particular attention in the recent decades due to the compelling demand for improving solar energy harvesting devices.^{1–7} In fact, these molecules, as well as chlorophyll (Chl) derivatives, act as reaction centers in both natural and artificial complex antenna systems.^{8–15} The key to understanding the initial steps of their photoexcited dynamics is the delicate interplay between the characteristic intense near-UV band (“*B*-band” or “Soret band,” around 400 nm) and the lower-energy visible band (“*Q*-band,” in the range of 500–600 nm),^{16} which are qualitatively understood in terms of Gouterman’s four-orbital model.^{17} The chemical functionalization of the free-base porphyrin (BP) generally preserves this excitation scheme, although it may affect the detailed shapes and positioning of the *Q* and *B* bands.^{3}

The dynamics between *B* and *Q* and within the *Q* band of BP and some derivatives has intensively been studied with both theoretical and experimental methods.^{18–23} The internal conversion times of BP in benzene solution were estimated^{18} to be 40 and 90 fs for the *B* → *Q*_{y} and *Q*_{y} → *Q*_{x} transition, respectively. In Ref. 19, it was reported that the charge transfer states (*CT*) appearing in diprotonated porphyrin may favor *B* → *Q*_{y} internal conversion through the indirect *B* → *CT* step. Time-resolved fluorescence experiments^{20} lead to the proposal of two different internal conversion pathways with different rates to explain the *B* band internal conversion in a tetra-phenyl-porphyrin, namely *B* → *Q*_{x} and *B* → *Q*_{y} → *Q*_{x}. In Ref. 22, linear response time-dependent density-functional theory (TDDFT)^{21} and on-the-fly fewest switches surface hopping (FSSH)^{24} were also employed to explain the relaxation process between the *B* and *Q* bands. It was shown therein that higher energy dark states (a band collectively called *N*) are involved in the *B* → *Q* internal conversion and that even *N* → *Q*_{x} population transfer is possible, even though less favorable than *N* → *Q*_{y}, provided that enough of excess energy is available. In Ref. 23, the FSSH approach was applied to describe non-radiative relaxation processes within the *Q*-bands of chlorophylls, showing the faster time crossing between the computed *Q*_{x} and *Q*_{y} population curves in the presence of the solvent phase as compared to the gas phase.

In the aforementioned studies, the internal dynamics following photo-excitation, relaxation times, and internal conversion pathways vary widely with respect to the ones of BP depending on the specific functionalization performed (for example, tetraphenylporphyrin,^{20} porphyrins bearing 0–4 meso-phenyl substituents^{2}). This fact renders each functionalized system unique and calls for theoretical characterization methods useful to find possible general trends.

Here, we focus on 5-ethoxycarbonyl-10-mesityl-15-benzyloxycarbonyl porphyrin^{25,26} (FP). The synthesis of this molecule involves placing a carboxylic acid group directly on one or more of the BP meso-carbon atoms.^{27} This allows for the construction of arrays in which the porphyrin macrocycles are close to each other and display an enhanced interaction with respect to, for example, the ones with hexa-phenylbenzene groups.^{28,29} We performed an in-depth analysis of the active normal modes, including an investigation of the potential energy surfaces (PESs) along their vibration trajectories, and compared the results obtained for both FP and BP. We show that excited-state normal-mode analysis, complemented by the calculation of per-mode reorganization energies (REs) and by a set of targeted scans along specific vibration modes, can unveil possible internal conversion pathways and point at specific regions of excited state PESs that can be crucial in non-adiabatic dynamics.

## II. METHODS

The analysis of normal modes on the excited state requires the calculation of per-mode reorganization energies (REs) and dimensionless Huang–Rhys (HR) factors, which provide a measure of the interaction strength between the electronic and vibrational states of the molecule. These quantities can be obtained within a displaced multi-mode harmonic oscillator model, which has successfully been applied in several studies.^{30–34} This approach, outlined below, is only rigorously valid as long as strong anharmonicities or Duschinsky effects^{35} can be neglected. We will also work in Condon excitation regime and neglect spin–orbit effects, confining ourselves to the singlet manifold.

*a*and

*b*, respectively, where we aim at determining the influence of the vibrational modes calculated on a selected state

*s*on the

*a*→

*b*transition. To define HRs and REs, we consider the shift of the potential energy surface (PES) Δ

*q*

_{μ}between

*a*and

*b*[see Fig. 1(a)], where the model represents the adiabatic harmonic potentials in the basis set of the normal modes

*μ*near the minimum of the PES of the selected state

*s*. Within the parallel harmonic approximation, the PES of a two-level system (such as in Fig. 1) for any value $q\mu *$ of the normal mode coordinate

*q*

_{μ}can be written in terms of displacements with respect to PES minima as

*a*(

*b*) at the minimum of its PES [see Fig. 1(a)], respectively, and

*M*

_{μ}and

*ω*

_{μ}are the modal mass and eigenfrequency, respectively.

*μ*[see Fig. 1(a)] as

*q*

_{μ}between the minima, the HR factor

*ξ*

_{μ}, and the RE

*E*

_{μ}for each mode

*μ*can then be written, respectively, as

^{36}

^{30–32,34,36}The spectral line shape

*I*(

*ω*) is defined in terms of the generation function

*G*(

*t*) as follows:

*T*is the temperature, and

*ξ*

_{μ}are the HR factors computed for the final state of the transition. The damping function is

^{30,37}In practice, since the high-frequency (“hard”) modes define the structure of the spectrum, while the low-frequency ones (“soft”) are responsible for the broadening, Γ can be calculated by splitting the normal modes into two groups and defining Γ as the average FWHM of the soft modes,

^{34,38}

Following the description of the adopted model, the actual procedure for computing HR factors [*ξ*_{μ}, Eq. (7)] and REs [*E*_{μ}, Eq. (8)] from quantum-mechanics (QM) calculations is outlined in Fig. 2 and fully described below. Specifically, the geometry of the initial state *a* is optimized, and it is then used to calculate the gradient of the total energy on the final state, *g*^{b}, at the coordinates corresponding to a vertical transition from *a*. The latter is obtained by computing the forces acting on each of the *N* atoms of the system in the state *b*. Then, the gradient *g*^{b} is projected onto the normal modes of the state *s*. It is worth noting that the choice of the state *s* for the calculation of the normal modes is usually dictated by the type of process under investigation, with *s* often chosen to coincide with *b*. Sometimes, another state could be chosen (e.g., the initial state or the ground state) as a cheaper approximation or, in case, a satisfying convergence on the excited state cannot be achieved, even though vibronic replica in the absorption spectra will likely be less precise in this case.^{39} Other possible choices of *s* are discussed in detail in Sec. III and in the supplementary material (Fig. S1).

*s*, we need to obtain its equilibrium geometry, for which the Hessian of the total energy has to be subsequently computed. From the diagonalization of the Hessian matrix, one can obtain the vibrational eigenfrequencies,

*ω*

_{μ}, and the reduced mass matrix, $M\u0302$, a (3

*N*− 6) × (3

*N*− 6) matrix whose diagonal elements are $1/M\mu $. Taking into account the relation between normal and Cartesian coordinates, i.e., $X=L\u0302M\u0302q$, where $L\u0302$ is the (3

*N*) × (3

*N*− 6) normalized transition matrix, we can then compute the gradient operator projections onto the normal modes,

*i*runs over the

*N*atoms of the system and

*j*runs over the three Cartesian components. Given $g\mu a,b$, one can then compute Δ

*q*

_{μ},

*ξ*

_{μ}, and

*E*

_{μ}for each normal mode

*μ*, according to Eqs. (6)–(8), respectively, and subsequently obtain the absorption spectrum with the inclusion of vibronic replicas, according to Eqs. (9)–(12).

In addition to providing spectra with vibronic effects, HR factors and REs are especially useful to deepen the analysis of the coupling with the nuclear degrees of freedom, as normal modes characterized by large HR factors and REs are more likely to play a role in the electronic transition of interest. Indeed, within the harmonic approximation, REs give an estimate of the energy variation along the excited-state PES [see Fig. 1(b)]. Once high-RE modes are identified, one can analyze the trajectories by moving the system along those specific modes. This kind of analysis is not meant to substitute explicitly dynamical methods,^{40} as the time variable does not appear. However, it can provide a simplified, intuitive picture of the adiabatic PES landscape for individual modes.

*μ*(Δ

*q*

_{μ}) in Cartesian coordinates (Δ

*X*

_{μ}) can be obtained back from the normal coordinates by using the vector of $L\u0302$ matrix along the

*μ*th mode. Therefore, the components of Δ

*X*

_{μ}vector are

*X*

_{init}, typically located at one PES bottom, deformed configurations following the

*μ*mode can then be computed as

*X*

_{fin}indicates the vector of the displaced coordinates along the normal mode, while

*X*

_{init}is the vector of the initial coordinates (the minimum of the initial state

*a*). At each desired

*X*

_{fin}, vertical excitation energies can be computed, possibly pinpointing “hot” regions where the PESs of different electronic states get critically close to each other or cross. An example of it is shown in Fig. 7.

## III. RESULTS AND DISCUSSION

Here, we consider, side by side, BP and FP. We apply the procedure detailed in Sec. II to characterize the PES landscape for modes actively taking part to both *B* → *Q* and *Q* → *Q* internal conversion processes. The QM calculations used to compute HR factors and REs are performed within the density-functional theory (DFT) and TDDFT frameworks by using the Gaussian 16 package.^{41} The hybrid range-corrected CAM-B3LYP^{42} functional, together with the 6-311(d,p) basis set, is used to determine both the equilibrium geometries and the gradients of the *a*, *b*, and *s* states. The effect of the solvent is included through the polarizable continuum model (PCM).^{43} The solvent used in the quoted experiments (see Table I) is tetrahydrofuran (THF) for BP^{16} and dichloromethane (DCM) for FP.^{25} Given the similar dielectric constant, THF is used for both BP and FP in TDDFT calculations to ease the comparison between the two systems. The optical transitions are characterized according to the natural transition orbital (NTO) analysis^{44} of the TDDFT transition density, as implemented in the Multiwfn package.^{45}

Band . | $BPexp16$ . | BP_{calc}
. | $FPexp25$ . | FP_{calc}
. | State . |
---|---|---|---|---|---|

Q_{x}(0–0) | 616 | 560 | 637 | 584 | S_{1} |

Q_{x}(0–1) | 561 | 582 | |||

Q_{y}(0–0) | 518 | 519 | 543 | 540 | S_{2} |

Q_{y}(0–1) | 487 | 506 | |||

B | 392 | 383/388 | 408 | 393/400 | S_{3}/S_{4} |

Band . | $BPexp16$ . | BP_{calc}
. | $FPexp25$ . | FP_{calc}
. | State . |
---|---|---|---|---|---|

Q_{x}(0–0) | 616 | 560 | 637 | 584 | S_{1} |

Q_{x}(0–1) | 561 | 582 | |||

Q_{y}(0–0) | 518 | 519 | 543 | 540 | S_{2} |

Q_{y}(0–1) | 487 | 506 | |||

B | 392 | 383/388 | 408 | 393/400 | S_{3}/S_{4} |

### A. Vibronic effects in absorption spectra

Figure 3 shows the calculated absorption spectra of both BP [black curve, panel (a)] and FP [blue curve, panel (b)]. The vertical lines indicate the zero-phonon excitations resulting from TDDFT simulations. Here, we follow the common convention and indicate the ground state as *S*_{0}. The (singlet) electronic excited states *S*_{1} and *S*_{2} belong to the Q band, while *S*_{3} and *S*_{4} belong to the B band. In particular, *S*_{1} and *S*_{2} correspond to the two non-degenerate *Q*_{x}(0–0) and *Q*_{y}(0–0) electronic excitations, arising from the lowered symmetry of BP with respect to metalloporphyrins^{46} (*D*_{2h} vs *D*_{4h}) due to the presence of NH protons. As we will see later in the discussion, this lowered symmetry not only affects the spectra but also plays an important role for the internal conversion dynamics.

A direct comparison between the electronic excitations for BP (a) and FP (b) shows that the excitation sequence remains the same upon functionalization, except for an overall redshift of the energies in the FP case, which are in quite good agreement with experimental data, as reported in Table I. The strength of the Q band features is larger in FP than in BP, in agreement with experimental data.^{16} In addition, the electronic excitations are accompanied by two phonon replicas each [*Q*_{x}(0–1), *Q*_{y}(0–1)], whose experimental values are also reported in Table I (see the discussion below).

Starting from the purely electronic spectra (Fig. 3, vertical bars), one can estimate the effect of molecular vibrations by computing the REs, as detailed in Sec. II. In the following, we focus on vibrational modes with high REs (hereafter called *active modes*), which are the ones contributing the most to the vibronic progression of the absorption spectrum and to the broadening of the peaks. Figure 4 displays the per-mode REs of BP (black bars) and FP (colored bars) for the transition from the ground to the different excited states, where the set of normal modes of *S*_{1} was used to compute the REs for the Q band transitions (from *S*_{0} to *S*_{1} and *S*_{2}) and the *S*_{3} set for the B band ones (from *S*_{0} to *S*_{3} and *S*_{4}). As can be noted by comparing the different panels of Fig. 4, the REs show an overall increase and a more spread distribution upon functionalization, irrespective of the chosen transition. In fact, BP (black bars) shows quite sparse and rather few active modes. On the contrary, FP shows a more spread “bath” of active vibrational modes, due to the further symmetry lowering caused by the presence of the functional groups attached to the core ring (for the atomic displacements along the active normal modes, see Fig. S2).

The absorption spectra computed by including vibronic effects are shown in Fig. 3. The vibration-induced broadening, Γ, is obtained as the mean FWHM of soft modes below 400 cm^{−1}. The difference between BP and FP in RE values and their distribution here yield different values for Γ, i.e., 514, 378, 458, and 380 cm^{−1} for the first four transitions of FP and, correspondingly, 389, 295, 279, and 352 cm^{−1} for BP. In the high-frequency-mode range, the different REs between BP and FP give rise, instead, to different vibronic progressions. Specifically, in the case of BP, a shoulder to the B band appears at ∼370 nm, which can be attributed to the 1777 cm^{−1} mode for the *S*_{0} → *S*_{3} transition and to the set of modes in the range 1400–1600 cm ^{−1} for the *S*_{0} → *S*_{4} transition. Vibronic replicas are, instead, almost negligible in the Q band due to the absence of modes with particularly high REs. Moving to FP, we find that the shoulder of the B band (∼378 nm) is noticeably more pronounced and originates from the contribution of the 2000 cm^{−1} mode for the *S*_{0} → *S*_{3} transition and from the set of modes at 1400–1600 cm^{−1} for the *S*_{0} → *S*_{4} transition. In addition, a vibronic peak at ∼510 nm arises in the Q band, mostly due to a high-RE mode at 1514 cm^{−1} for the *S*_{0} → *S*_{2} transition. The position of the *S*_{2} vibronic peak is close to the experimental *Q*_{y}(0–1) (Table I), while it is impossible to clearly define a vibronic peak corresponding to *Q*_{x}(0–1) due to slightly blueshifted electronic peak *S*_{1} [*Q*_{x}(0–0)]. Indeed, as already shown for BP,^{47,48} it is known that the inclusion of Herzberg–Teller effect is important to better reproduce absorption line shapes, which is, however, beyond the scope of this work.

### B. PES along active normal modes

In addition to correcting the absorption spectra for vibronic effects, the per-mode REs are especially useful to understand the relaxation and internal conversion pathways. We have thus further examined the per-mode REs for different excited-state transitions, both within the Q band and between B and Q bands [see Figs. 5(a) and 5(b), respectively]. Here, REs are calculated on the set of *S*_{1} modes. Again, a remarkable difference is found between BP (in black), with few sparse and weakly active modes, and FP (in colors), with several active modes, especially in the range above 1200 cm^{−1}. By focusing on the most active vibrations (highlighted with arrows in Fig. 5) and by inspecting the corresponding atomic displacements (see Fig. 6), we find that these modes mostly differ in character for BP and FP, despite being all in-plane modes. In particular, the 1270 cm^{−1} mode of FP involves large motion of the double carbon bonds of the ethoxy-carbonyl and benzyloxy-carbonyl connectors. Moreover, both the modes at 1270 cm^{−1} and 1370 cm^{−1} show asymmetric rocking of the N–H groups, not seen in the case of BP. Only the vibration at about 1500 cm^{−1} has a similar pattern in the two molecules, although the amplitudes are less symmetric for FP. Notably, the latter mode has high RE for transitions both within the Q band and between Q and B bands, for both molecules. An in-depth analysis of the PES along this and other active modes by scanning vibrational trajectories can thus provide valuable information, detailed below.

#### 1. Single-mode analysis

We have, so far, separately determined the electronic structure and the active vibrational modes on the excited states of BP and FP. Now, we can merge this information by reconstructing the PES of the systems along the selected active modes. Figure 7 shows the energy-level scheme (a) and the Kohn–Sham molecular orbitals (MOs) of the optimized *S*_{1} state for BP (b) and FP (c). The HOMO and HOMO-1 orbitals look alike in the two molecules, except that the symmetries of BP HOMO-1 and HOMO (*A*_{u} and *B*_{1u}, respectively) are exchanged in FP, and the appearance of some density localized on the O atoms of the carboxylic acid group in FP. The LUMO and LUMO+1 MOs are, instead, remarkably affected by the connectors, which lower the symmetry and allow for a different mixing of the states. Moreover, by inspecting a few more states that are involved in higher-energy excitations [see Figs. 7(d) and 7(e) and discussion below], we find that HOMO-2 and HOMO-3 in BP [Fig. 7(b)] correspond to HOMO-4 and HOMO-5 in FP [Fig. 7(c)]; HOMO-2 and HOMO-3 in FP [Fig. 7(c)] are, instead, completely localized on the mesityl group.

Starting from the above analysis and using Eqs. (15) and (16), we define displaced geometries and compute the PESs along the highest RE modes at 1541 cm^{−1} in BP and at 1514 cm^{−1} in FP [see Fig. 5(a)], which were found to have similar characters in the two molecules (see Fig. 6 and the discussion above). Figure 7 shows a cut of the PESs along the two selected modes for BP [panel (d)] and FP [panel (e)]. The PESs, calculated at the dotted points, are interpolated according to the harmonic approximation; their colors refer to the color code of the group of occupied MOs involved in the transitions [Figs. 7(a)–7(c)]. This analysis allows one to understand whether the transition is dominated by Gouterman’s “dynamics” (red) or other orbitals are involved (orange and green). Notably, there is no influence of the orbitals localized on mesityl group (green ones) on the B band, whereas they contribute to higher excited states (see Fig. S3). Meanwhile, even though the scan only represents a specific section of the actual multidimensional space, it clearly shows that the *S*_{5} and *S*_{6} states (we label them collectively as *N*, as in Ref. 22) are crossing the B band. This is consistent with the earlier suggested mechanism where upper energy levels favor *B* → *Q*_{y} internal conversion through the indirect *B* → *N* step.^{19,22}

In addition to the crossing between *B* and *N* states, by looking at the reconstructed PESs along the selected mode, we notice the existence of a point toward which the gap between the *S*_{1} and *S*_{2} PESs tends to vanish. While this happens in both molecules, the energetics is rather different in the two cases. In fact, the excess energy of the crossing point $(q\u0304)$ with respect to the *S*_{2} minimum $\Delta Eexc=ES2(q\u0304)\u2212ES2(q0)$ is 1.1 eV in the case of BP and 0.31 eV in FP. The same occurs with respect to the *S*_{1} minimum. We conclude that the crossing point between the two Q-band states, driven by the selected high-frequency mode (1541 cm^{−1} in BP and 1514 cm^{−1} in FP), is much more easily accessible in FP than in BP. The fact that similar modes appear in both BP and FP, but with greatly enhanced REs in the latter, suggests the existence of a measurable effect in the intra-band dynamics, such as a much faster internal conversion time.

Let us now examine more closely the Q band of FP [Fig. 8(a)], by zooming in the region of the PES, where the *S*_{1} to *S*_{2} gap vanishes [panels (b) and (d)]. To identify the nature of the states on each surface around the zero-gap point, we considered three points along the trajectory centered around Δ*q* ≈ −15.7, namely A, B, and C [Fig. 8(a)], where the numbers (1 or 2) indicate the order of the excited state, i.e., *S*_{1} or *S*_{2}. For these selected q coordinates, we have analyzed both the oscillator strength [panel (b)] and the transition dipole moments [(TDM), panel (d)] of *S*_{1} and *S*_{2} along the symmetry directions of BP [x and y axes, as defined in panel (c) and corresponds to the states of *Q*_{x} and *Q*_{y} bands in Table I]. The oscillator strength, represented by the line thickness in panel (b), decreases but remains non-zero moving from *A*1 to *C*2, while it remains closer to zero moving from *A*2 to *C*1. The TDM display a similar trend, with x (black arrows) prevailing component (*C*2 → *B*1 → *A*1) or y (blue arrows) one (*C*1 → *B*2 → *A*2), again suggesting a crossing of the states.

For the same selected states, we also computed the NTOs in order to analyze their composition and character (see details in Sec. II). Specifically, we used the NTOs of *S*_{1} at point A (*A*1)—here the Gouterman MOs—as the basis set for our analysis, which is reported in Fig. 8(e). The gray scale indicates the composition of each selected state with respect to *A*1, partitioned onto its NTOs. A diagonal pattern indicates that the nature/character of the selected state is the same as *A*1; the presence of off-diagonal elements indicates, instead, that the nature of the state is different from that of the reference one. For instance, projecting the *A*2 state on *A*1 shows the exchange of HOMO and HOMO-1 orbitals, in accordance with the nature of *Q*_{x} and *Q*_{y} as described by the Gouterman model.^{17} As for oscillator strength and TDM trends, also the pattern found by the NTO analysis points to a crossing of the *S*_{1} and *S*_{2} states.

All of these analyses allow us to closely follow the character of the states around the zero-gap point of *S*_{1} and *S*_{2}. This points to an actual exchange of the characters, instead of a repulsion of the PESs. As the adiabatic approximation holds only far from this point, it is, therefore, plausible to assume that the two surfaces will actually give rise to a conical intersection or, at least, will become non-adiabatically coupled in the neighborhood of the crossing point. However, the investigation of the detailed geometry of the two surfaces and the classification of the intersection require dedicated methods, beyond the domain of the approximations adopted here, and are left for future investigation.

#### 2. Multi-modal analysis

In order to better understand the role of vibrations in the relaxation dynamics, we have explored the PES along additional modes (see Fig. 9), which could influence the *S*_{2} → *S*_{1} internal conversion by acting cooperatively with the highest-RE mode analyzed previously. For FP, we explore the 2D space defined by the highest-RE mode (*main*, 1514 cm^{−1}, RE = 75.08 cm^{−1}) with three other modes in the same high-frequency region, having *low* [1348 cm^{−1}, RE = 0.95 cm^{−1}, panel (a)], *medium* [1234 cm^{−1}, RE = 34.08 cm^{−1}, panel (b)], and *high* [1270 cm^{−1}, RE = 52.29 cm^{−1}, panel (c)] REs, respectively; for BP, we combine the highest-RE mode (1541 cm^{−1}, RE = 5.04 cm^{−1}) with the next-highest-RE one [1366 cm^{−1}, RE = 1.79 cm^{−1}, panel (d)]. The color maps reported in Fig. 9 display the absolute value of the energy difference |Δ*E*| between *S*_{2} and *S*_{1} in the 2D manifold defined by the two selected modes. This analysis provides valuable information, not only on the existence, location, and shape of critical points/regions where the system can display a strong non-adiabatic coupling (brown to black areas) but also on the energetic accessibility of these points, e.g., Δ*E*_{exc} between the *S*_{2} minimum and the touching point $(q\u0304)$.

By looking at the different 2D maps computed for FP, we can notice that the weakly active mode at 1348 cm^{−1} does not have any cooperative effect. In fact, the touching region is almost parallel to the horizontal axis in the plot, that is, Δ*E* remains nearly the same by moving along the normal coordinate of this *low*-RE mode [Fig. 9(a)]. On the contrary, modes with higher RE values [panels (b) and (c)] can significantly modify the local landscape, leading to smaller excess energy from the *S*_{2} minimum. The *medium*-RE mode [panel (b)] has Δ*E*_{exc} = 0.19 eV at the crossing point on the PES [$q\u0304=(\u22126.0,\u221212.7)$, marked with a white cross]; the *high*-RE mode [panel (c)] leads to Δ*E*_{exc} = 0.16 eV $[q\u0304=(\u22129.8,\u221210.2)]$.

In conclusion, the comparison of the maps obtained for FP and BP [Figs. 9(c) and 9(d)] clearly shows that the crossing region is much more accessible for FP than for BP, which confirms the possibility of a faster relaxation in FP, as anticipated from the single-mode analysis. In fact, for BP, we have found a barrier of 0.92 eV [panel (d), $q\u0304=(14.0,\u221218.0)$], whereas the excess energy for FP given by the cooperative effect of the two highest-RE modes is four times lower, i.e., Δ*E*_{exc} = 0.16 eV.

## IV. CONCLUSIONS

We have employed excited-state normal mode analysis to explore the PESs of BP and FP. This involved taking the following steps: (1) defining active modes (i.e., finding the particular set of modes of interest and selecting the modes with the highest RE values); (2) building transition energy scans along the active modes of interest; and (3) analyzing the states near the critical regions using the trends in changing oscillator strengths and transition dipole moments (values and directions) and by comparing natural transition orbitals between the states of different structures along the normal mode scans.

All of these analyses point at a crossing between the PESs of *Q*_{x} and *Q*_{y} states and suggest that the considered functionalization of the porphyrin may substantially enhance the internal conversion within the *Q* band. Moreover, the study of the 2D PES along the active modes demonstrates that the FP has a much higher probability than the BP to reach the crossing point between the *Q*_{x} and *Q*_{y} states, upon *Q* band excitation. The barrier between *Q*_{y} minimum and the crossing point is just 0.16 eV for FP, whereas it is 0.92 eV for BP.

We must bear in mind that the methodology followed here is bound to several constraints, as it is only rigorously valid in the adiabatic regime, within the assumption of harmonic PES for each of the considered vibrational modes and excluding vibrational mixing. However, the method can be profitably exploited as a tool to quickly scan the excited-state PES for interesting points—following the lead of the most active modes—at a reasonable computational cost. Therefore, it can be used as a convenient initial step for exploring the effects of functionalization on the internal dynamics of porphyrins and other molecules. In the specific case examined here, we theoretically explain how a specific functionalization may have a huge impact on the internal conversion within the *Q* band. We identify a particular vibrational mode responsible for driving the system into a conversion sweet-point, which opens the road, on the one hand, to a more advanced investigation of the dynamics at or close to the critical region and, on the other hand, to a more systematic study of different functionalization schemes.

## SUPPLEMENTARY MATERIAL

Possible choices for the set of vibrational modes for HR and RE calculations, as well as atomic displacements for the modes with the highest RE values for the transitions in Fig. 4 and the orbital composition of the excited states along the trajectories shown in Figs. 7(d) and 7(e), can be found in the supplementary material.

## ACKNOWLEDGMENTS

This work was supported by the Italian Ministry of University and Research, within the program PRIN 2017, Grant No. 201795SBA3—HARVEST. The computational time on the Marconi100 machine at CINECA was provided by the Italian ISCRA program.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Pavel Rukin**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Deborah Prezzi**: Conceptualization (equal); Funding acquisition (lead); Project administration (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal). **Carlo Andrea Rozzi**: Conceptualization (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## REFERENCES

*meso*-phenyl substituents: New insights into seemingly well understood tetrapyrroles

*π*-expanded chlorophyll derivatives to construct light-harvesting antenna models

*Chl*

_{2A}

*P*

_{A}

*P*

_{B}

*Chl*

_{2B}chlorophyll exciplex in photosystem I

*Time-Dependent Density-Functional Theory: Concepts and Applications*

*π–π*interchromophore interactions

^{4}T

_{1}→

^{6}A

_{1}emission bands of Mn

^{2+}impurity in Zn

_{2}SiO

_{4}: A first-principle methodology

^{2+}ions in ZnCo

_{3}and CaCo

_{3}matrices by means of quantum chemistry

*ab initio*calculations

_{1}and S

_{2}states of pyridine

*Q*

_{x}band of porphyrin as a case study

_{x}band in free-base porphin