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.

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 estimated18 to be 40 and 90 fs for the BQy and QyQx transition, respectively. In Ref. 19, it was reported that the charge transfer states (CT) appearing in diprotonated porphyrin may favor BQy internal conversion through the indirect BCT step. Time-resolved fluorescence experiments20 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 BQx and BQyQx. 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 BQ internal conversion and that even NQx population transfer is possible, even though less favorable than NQy, 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 Qx and Qy 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 substituents2). 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 porphyrin25,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.

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 effects35 can be neglected. We will also work in Condon excitation regime and neglect spin–orbit effects, confining ourselves to the singlet manifold.

We are normally concerned about transitions (could be either optical absorption or internal conversion) occurring between an initial state and a final electronic state, hereafter labeled a and b, respectively, where we aim at determining the influence of the vibrational modes calculated on a selected state s on the ab 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μ* of the normal mode coordinate qμ can be written in terms of displacements with respect to PES minima as
(1)
where Ea0 (Eb0) and qμa0 (qμb0) are the energies and the normal coordinates of the initial (final) state a (b) at the minimum of its PES [see Fig. 1(a)], respectively, and Mμ and ωμ are the modal mass and eigenfrequency, respectively.
FIG. 1.

Definition of the shift Δqμ along the mode μ in normal coordinates between an initial state a and a final state b. (a) General case with an arbitrary chosen value of the normal coordinate qμ*; (b) special common case in which qμ* is chosen at the minimum of the PES of state a. Eμa and Eμb are the PESs of states a and b along the mode μ, respectively; qμa0, qμb0 and Eμa0, Eμb0 are the coordinates and energies of the minima of the initial and final states, respectively; Eμan and Eμbm are the energies of the vibrationally excited states, respectively; gμa,gμb are the gradients of a and b PESs at the qμ* normal coordinate, respectively.

FIG. 1.

Definition of the shift Δqμ along the mode μ in normal coordinates between an initial state a and a final state b. (a) General case with an arbitrary chosen value of the normal coordinate qμ*; (b) special common case in which qμ* is chosen at the minimum of the PES of state a. Eμa and Eμb are the PESs of states a and b along the mode μ, respectively; qμa0, qμb0 and Eμa0, Eμb0 are the coordinates and energies of the minima of the initial and final states, respectively; Eμan and Eμbm are the energies of the vibrationally excited states, respectively; gμa,gμb are the gradients of a and b PESs at the qμ* normal coordinate, respectively.

Close modal
We obtain the projections of the gradients of the PES along the normal mode μ [see Fig. 1(a)] as
(2)
The normal coordinate displacement Δqμ between the minima, the HR factor ξμ, and the RE Eμ for each mode μ can then be written, respectively, as36 
(3)
(4)
(5)
Within the harmonic approximation, the difference of gradients gμbgμa does not depend on the initial point qμ*. As such, qμ* can be chosen in the most convenient way, e.g., such as to minimize computing time. If qμ* is taken as the minimum of the initial state qμa0 [see Fig. 1(b)], the coordinate differences, HR factors, and REs assume the simplified form,
(6)
(7)
(8)
Furthermore, from the knowledge of HR factors [Eq. (7)], one can obtain the absorption spectrum with the inclusion of vibronic replicas (within the Franck–Condon approximation) by using the generating function approach.30–32,34,36 The spectral line shape I(ω) is defined in terms of the generation function G(t) as follows:
(9)
(10)
where Ea0Eb0 is the purely electronic (zero-phonon) transition energy, T is the temperature, and ξμ are the HR factors computed for the final state of the transition. The damping function is
(11)
where Γ is the homogeneous linewidth.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
(12)
where
(13)

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, gb, 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 gb 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).

FIG. 2.

Flow chart of the procedure used to compute HR factors and REs, as described in the text.

FIG. 2.

Flow chart of the procedure used to compute HR factors and REs, as described in the text.

Close modal
For a given choice of 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̂, a (3N − 6) × (3N − 6) matrix whose diagonal elements are 1/Mμ. Taking into account the relation between normal and Cartesian coordinates, i.e., X=L̂M̂q, where L̂ is the (3N) × (3N − 6) normalized transition matrix, we can then compute the gradient operator projections onto the normal modes,
(14)
where i runs over the N atoms of the system and j runs over the three Cartesian components. Given gμ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.

The displacement along a mode μqμ) in Cartesian coordinates (ΔXμ) can be obtained back from the normal coordinates by using the vector of L̂ matrix along the μth mode. Therefore, the components of ΔXμ vector are
(15)
Starting from the initial configuration Xinit, typically located at one PES bottom, deformed configurations following the μ mode can then be computed as
(16)
where Xfin indicates the vector of the displaced coordinates along the normal mode, while Xinit is the vector of the initial coordinates (the minimum of the initial state a). At each desired Xfin, 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.

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 BQ and QQ 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-B3LYP42 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 BP16 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) analysis44 of the TDDFT transition density, as implemented in the Multiwfn package.45 

TABLE I.

Comparison between experimental data16,25 and calculated vertical transitions (in nm).

BandBPexp16BPcalcFPexp25FPcalcState
Qx(0–0) 616 560 637 584 S1 
Qx(0–1) 561  582   
Qy(0–0) 518 519 543 540 S2 
Qy(0–1) 487  506   
392 383/388 408 393/400 S3/S4 
BandBPexp16BPcalcFPexp25FPcalcState
Qx(0–0) 616 560 637 584 S1 
Qx(0–1) 561  582   
Qy(0–0) 518 519 543 540 S2 
Qy(0–1) 487  506   
392 383/388 408 393/400 S3/S4 

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 S0. The (singlet) electronic excited states S1 and S2 belong to the Q band, while S3 and S4 belong to the B band. In particular, S1 and S2 correspond to the two non-degenerate Qx(0–0) and Qy(0–0) electronic excitations, arising from the lowered symmetry of BP with respect to metalloporphyrins46 (D2h vs D4h) 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.

FIG. 3.

Calculated linear absorption spectra of BP in black [panel (a)] and FP in blue [panel (b)], including both vibronic effects and broadening. The vertical bars indicate the vertical transition energies normalized to the maximum of the oscillator strength. Panels (c) and (d) show the ground-state structures of BP (black) and FP (blue).

FIG. 3.

Calculated linear absorption spectra of BP in black [panel (a)] and FP in blue [panel (b)], including both vibronic effects and broadening. The vertical bars indicate the vertical transition energies normalized to the maximum of the oscillator strength. Panels (c) and (d) show the ground-state structures of BP (black) and FP (blue).

Close modal

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 [Qx(0–1), Qy(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 S1 was used to compute the REs for the Q band transitions (from S0 to S1 and S2) and the S3 set for the B band ones (from S0 to S3 and S4). 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).

FIG. 4.

Calculated per-mode reorganization energies (REs) based on transitions from the ground state (S0) to (a) S1 and (b) S2 (Q band) and (c) S3 and (d) S4 (B band) excited states. REs of FP modes are in colors (red for S1, green for S2, blue for S3, and magenta for S4); BP REs are in black. The PES depicted in dark red in the schematic indicates the state chosen to compute the basis set of the normal modes for RE calculations.

FIG. 4.

Calculated per-mode reorganization energies (REs) based on transitions from the ground state (S0) to (a) S1 and (b) S2 (Q band) and (c) S3 and (d) S4 (B band) excited states. REs of FP modes are in colors (red for S1, green for S2, blue for S3, and magenta for S4); BP REs are in black. The PES depicted in dark red in the schematic indicates the state chosen to compute the basis set of the normal modes for RE calculations.

Close modal

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 S0S3 transition and to the set of modes in the range 1400–1600 cm −1 for the S0S4 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 S0S3 transition and from the set of modes at 1400–1600 cm−1 for the S0S4 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 S0S2 transition. The position of the S2 vibronic peak is close to the experimental Qy(0–1) (Table I), while it is impossible to clearly define a vibronic peak corresponding to Qx(0–1) due to slightly blueshifted electronic peak S1 [Qx(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.

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 S1 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.

FIG. 5.

Calculated per-mode REs for (a) transitions from S1 to S2 (intra-Q band) and (b) from S2 to S3 (from Q to B band) excited states. REs of FP modes are in colors (teal for S1 to S2 and olive for S2 to S3; BP REs are in black. The PES depicted in dark red in the schematic indicates the state chosen to compute the basis set of the normal modes for RE calculations. The arrows and values highlight the most active modes.

FIG. 5.

Calculated per-mode REs for (a) transitions from S1 to S2 (intra-Q band) and (b) from S2 to S3 (from Q to B band) excited states. REs of FP modes are in colors (teal for S1 to S2 and olive for S2 to S3; BP REs are in black. The PES depicted in dark red in the schematic indicates the state chosen to compute the basis set of the normal modes for RE calculations. The arrows and values highlight the most active modes.

Close modal
FIG. 6.

Atomic displacements for the normal modes with higher REs in the S1S2 transition for (a) BP and (b) FP.

FIG. 6.

Atomic displacements for the normal modes with higher REs in the S1S2 transition for (a) BP and (b) FP.

Close modal

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 S1 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 (Au and B1u, 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.

FIG. 7.

(a) Energy-level scheme for BP and FP. (b) and (c) Selected molecular orbitals (MOs), contributing to the lowest excited state transitions, for BP and FP, respectively. (d) and (e) Trajectory scan along the modes at 1541 cm−1 in BP and 1514 cm−1 in FP, respectively. The q = 0 coordinate corresponds to the optimized geometry in S1, for which energy levels (a) and MOs [(b) and (c)] are shown. The excited-state PESs are colored according to the color code of the MOs [(b) and (c)], contributing the most to the transitions. The ground-state PES is in gray. The dots represent the actual vertical energy calculations; the solid lines are obtained by interpolation, according to the harmonic approximation.

FIG. 7.

(a) Energy-level scheme for BP and FP. (b) and (c) Selected molecular orbitals (MOs), contributing to the lowest excited state transitions, for BP and FP, respectively. (d) and (e) Trajectory scan along the modes at 1541 cm−1 in BP and 1514 cm−1 in FP, respectively. The q = 0 coordinate corresponds to the optimized geometry in S1, for which energy levels (a) and MOs [(b) and (c)] are shown. The excited-state PESs are colored according to the color code of the MOs [(b) and (c)], contributing the most to the transitions. The ground-state PES is in gray. The dots represent the actual vertical energy calculations; the solid lines are obtained by interpolation, according to the harmonic approximation.

Close modal

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 S5 and S6 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 BQy internal conversion through the indirect BN 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 S1 and S2 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̄) with respect to the S2 minimum ΔEexc=ES2(q̄)ES2(q0) is 1.1 eV in the case of BP and 0.31 eV in FP. The same occurs with respect to the S1 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 S1 to S2 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., S1 or S2. For these selected q coordinates, we have analyzed both the oscillator strength [panel (b)] and the transition dipole moments [(TDM), panel (d)] of S1 and S2 along the symmetry directions of BP [x and y axes, as defined in panel (c) and corresponds to the states of Qx and Qy bands in Table I]. The oscillator strength, represented by the line thickness in panel (b), decreases but remains non-zero moving from A1 to C2, while it remains closer to zero moving from A2 to C1. The TDM display a similar trend, with x (black arrows) prevailing component (C2 → B1 → A1) or y (blue arrows) one (C1 → B2 → A2), again suggesting a crossing of the states.

FIG. 8.

(a) PES of the Q band of FP along the mode at 1514 cm−1 [S0 (gray), S1 and S2 (red)]. The solid dots represent the actual vertical energy calculations; the solid lines are obtained by interpolation, according to the harmonic approximation. A1, A2, B1, B2, C1, and C2 are selected transitions near the possible crossing, where the numbers indicate the excited state order, i.e., S1 or S2. (b) Zoomed-in view of the possible crossing displayed in (a). Oscillator strengths are represented by the line thickness, while the numerical values are indicated for the selected transitions. (d) Scan of the transition dipole moment (TDM) direction, with respect to the axes indicated in panel (c) (X—black, Y—blue), chosen as the symmetry directions of BP. The length of the arrows represents the TDM value. Energy values in panels (b) and (d) are shown with respect to S0 at each normal coordinate (ESiES0). (e) NTO composition of the selected transitions (represented on the y axis) with respect to a reference state, here chosen to be A1 (x axis). The intensity of the color indicates the weights of A1 NTOs with respect to the Gouterman MOs.

FIG. 8.

(a) PES of the Q band of FP along the mode at 1514 cm−1 [S0 (gray), S1 and S2 (red)]. The solid dots represent the actual vertical energy calculations; the solid lines are obtained by interpolation, according to the harmonic approximation. A1, A2, B1, B2, C1, and C2 are selected transitions near the possible crossing, where the numbers indicate the excited state order, i.e., S1 or S2. (b) Zoomed-in view of the possible crossing displayed in (a). Oscillator strengths are represented by the line thickness, while the numerical values are indicated for the selected transitions. (d) Scan of the transition dipole moment (TDM) direction, with respect to the axes indicated in panel (c) (X—black, Y—blue), chosen as the symmetry directions of BP. The length of the arrows represents the TDM value. Energy values in panels (b) and (d) are shown with respect to S0 at each normal coordinate (ESiES0). (e) NTO composition of the selected transitions (represented on the y axis) with respect to a reference state, here chosen to be A1 (x axis). The intensity of the color indicates the weights of A1 NTOs with respect to the Gouterman MOs.

Close modal

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 S1 at point A (A1)—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 A1, partitioned onto its NTOs. A diagonal pattern indicates that the nature/character of the selected state is the same as A1; 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 A2 state on A1 shows the exchange of HOMO and HOMO-1 orbitals, in accordance with the nature of Qx and Qy 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 S1 and S2 states.

All of these analyses allow us to closely follow the character of the states around the zero-gap point of S1 and S2. 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 S2S1 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 S2 and S1 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., ΔEexc between the S2 minimum and the touching point (q̄).

FIG. 9.

Absolute value of the energy difference |ΔE| between S2 and S1 in FP [(a)–(c)] and BP (d), obtained by exploring the PES along two selected modes, i.e., the one with the highest RE (main, y axis) and a second mode (x axis) chosen to have low (a), medium (b), or high [(c) and (d)] RE. The axes show the displacement along the modes in normal coordinates, where q = 0 corresponds to the S1 optimized geometry. The white contour lines and the values (in eV) indicate the energy of S2 PES relatively to the minimum of the S0 PES. The white circle and the related value (in eV) indicate the minimum of S2 along the selected modes. The white cross indicates the point with the lowest |ΔE| value.

FIG. 9.

Absolute value of the energy difference |ΔE| between S2 and S1 in FP [(a)–(c)] and BP (d), obtained by exploring the PES along two selected modes, i.e., the one with the highest RE (main, y axis) and a second mode (x axis) chosen to have low (a), medium (b), or high [(c) and (d)] RE. The axes show the displacement along the modes in normal coordinates, where q = 0 corresponds to the S1 optimized geometry. The white contour lines and the values (in eV) indicate the energy of S2 PES relatively to the minimum of the S0 PES. The white circle and the related value (in eV) indicate the minimum of S2 along the selected modes. The white cross indicates the point with the lowest |ΔE| value.

Close modal

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 S2 minimum. The medium-RE mode [panel (b)] has ΔEexc = 0.19 eV at the crossing point on the PES [q̄=(6.0,12.7), marked with a white cross]; the high-RE mode [panel (c)] leads to ΔEexc = 0.16 eV [q̄=(9.8,10.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̄=(14.0,18.0)], whereas the excess energy for FP given by the cooperative effect of the two highest-RE modes is four times lower, i.e., ΔEexc = 0.16 eV.

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 Qx and Qy 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 Qx and Qy states, upon Q band excitation. The barrier between Qy 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.

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.

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.

The authors have no conflicts to disclose.

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).

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

1.
S.
Takagi
,
M.
Eguchi
,
D. A.
Tryk
, and
H.
Inoue
, “
Porphyrin photochemistry in inorganic/organic hybrid materials: Clays, layered semiconductors, nanotubes, and mesoporous materials
,”
J. Photochem. Photobiol., C
7
,
104
126
(
2006
).
2.
A. K.
Mandal
,
M.
Taniguchi
,
J. R.
Diers
,
D. M.
Niedzwiedzki
,
C.
Kirmaier
,
J. S.
Lindsey
,
D. F.
Bocian
, and
D.
Holten
, “
Photophysical properties and electronic structure of porphyrins bearing zero to four meso-phenyl substituents: New insights into seemingly well understood tetrapyrroles
,”
J. Phys. Chem. A
120
,
9719
9731
(
2016
).
3.
S.
Hiroto
,
Y.
Miyake
, and
H.
Shinokubo
, “
Synthesis and functionalization of porphyrins through organometallic methodologies
,”
Chem. Rev.
117
,
2910
3043
(
2017
).
4.
R.
Paolesse
,
S.
Nardis
,
D.
Monti
,
M.
Stefanelli
, and
C.
Di Natale
, “
Porphyrinoids for chemical sensor applications
,”
Chem. Rev.
117
,
2517
2583
(
2017
).
5.
M. O.
Senge
,
N. N.
Sergeeva
, and
K. J.
Hale
, “
Classic highlights in porphyrin and porphyrinoid total synthesis and biosynthesis
,”
Chem. Soc. Rev.
50
,
4730
4789
(
2021
).
6.
M. K.
Panda
,
K.
Ladomenou
, and
A. G.
Coutsolelos
, “
Porphyrins in bio-inspired transformations: Light-harvesting to solar cell
,”
Coord. Chem. Rev.
256
,
2601
2627
(
2012
).
7.
C.
Biswas
,
S. G.
Palivela
,
L.
Giribabu
,
V. R.
Soma
, and
S. S. K.
Raavi
, “
Femtosecond excited-state dynamics and ultrafast nonlinear optical investigations of ethynylthiophene functionalized porphyrin
,”
Opt. Mater.
127
,
112232
(
2022
).
8.
J. G.
Woller
,
J. K.
Hannestad
, and
B.
Albinsson
, “
Self-assembled nanoscale DNA–porphyrin complex for artificial light harvesting
,”
J. Am. Chem. Soc.
135
,
2759
2768
(
2013
).
9.
J.
Otsuki
, “
Supramolecular approach towards light-harvesting materials based on porphyrins and chlorophylls
,”
J. Mater. Chem. A
6
,
6710
6753
(
2018
).
10.
S.
Matsubara
and
H.
Tamiaki
, “
Synthesis and self-aggregation of π-expanded chlorophyll derivatives to construct light-harvesting antenna models
,”
J. Org. Chem.
83
,
4355
4364
(
2018
).
11.
W.
Auwärter
,
D.
Écija
,
F.
Klappenberger
, and
J. V.
Barth
, “
Porphyrins at interfaces
,”
Nat. Chem.
7
,
105
120
(
2015
).
12.
M. J.
Llansola-Portoles
,
D.
Gust
,
T. A.
Moore
, and
A. L.
Moore
, “
Artificial photosynthetic antennas and reaction centers
,”
C. R. Chim.
20
,
296
313
(
2017
).
13.
M.
Gorka
,
P.
Charles
,
V.
Kalendra
,
A.
Baldansuren
,
K.
Lakshmi
, and
J. H.
Golbeck
, “
A dimeric chlorophyll electron acceptor differentiates type I from type II photosynthetic reaction centers
,”
iScience
24
,
102719
(
2021
).
14.
V.
Mascoli
,
V.
Novoderezhkin
,
N.
Liguori
,
P.
Xu
, and
R.
Croce
, “
Design principles of solar light harvesting in plants: Functional architecture of the monomeric antenna CP29
,”
Biochim. Biophys. Acta, Bioenerg.
1861
,
148156
(
2020
).
15.
D. A.
Cherepanov
,
I. V.
Shelaev
,
F. E.
Gostev
,
A.
Petrova
,
A. V.
Aybush
,
V. A.
Nadtochenko
,
W.
Xu
,
J. H.
Golbeck
, and
A. Y.
Semenov
, “
Primary charge separation within the structurally symmetric tetrameric Chl2APAPBChl2B chlorophyll exciplex in photosystem I
,”
J. Photochem. Photobiol., B
217
,
112154
(
2021
).
16.
J.
Braun
,
C.
Hasenfratz
,
R.
Schwesinger
, and
H.-H.
Limbach
, “
Free acid porphyrin and its conjugated monoanion
,”
Angew Chem. Int. Ed. Engl.
33
,
2215
2217
(
1994
).
17.
M.
Gouterman
, “
Spectra of porphyrins
,”
J. Mol. Spectrosc.
6
,
138
163
(
1961
).
18.
S.
Akimoto
,
T.
Yamazaki
,
I.
Yamazaki
, and
A.
Osuka
, “
Excitation relaxation of zinc and free-base porphyrin probed by femtosecond fluorescence spectroscopy
,”
Chem. Phys. Lett.
309
,
177
182
(
1999
).
19.
A.
Marcelli
,
P.
Foggi
,
L.
Moroni
,
C.
Gellini
, and
P. R.
Salvi
, “
Excited-state absorption and ultrafast relaxation dynamics of porphyrin, diprotonated porphyrin, and tetraoxaporphyrin dication
,”
J. Phys. Chem. A
112
,
1864
1872
(
2008
).
20.
S. Y.
Kim
and
T.
Joo
, “
Coherent nuclear wave packets in Q states by ultrafast internal conversions in free base tetraphenylporphyrin
,”
J. Phys. Chem. Lett.
6
,
2993
2998
(
2015
).
21.
C. A.
Ullrich
,
Time-Dependent Density-Functional Theory: Concepts and Applications
(
Oxford University Press
,
Oxford
,
2011
).
22.
K.
Falahati
,
C.
Hamerla
,
M.
Huix-Rotllant
, and
I.
Burghardt
, “
Ultrafast photochemistry of free-base porphyrin: A theoretical investigation of B → Q internal conversion mediated by dark states
,”
Phys. Chem. Chem. Phys.
20
,
12483
12492
(
2018
).
23.
M.
Fortino
,
E.
Collini
,
J.
Bloino
, and
A.
Pedone
, “
Unraveling the internal conversion process within the Q-bands of a chlorophyll-like-system through surface-hopping molecular dynamics simulations
,”
J. Chem. Phys.
154
,
094110
(
2021
).
24.
J. E.
Subotnik
,
W.
Ouyang
, and
B. R.
Landry
, “
Can we derive Tully’s surface-hopping algorithm from the semiclassical quantum Liouville equation? Almost, but only with decoherence
,”
J. Chem. Phys.
139
,
214107
(
2013
).
25.
Y.
Terazono
,
G.
Kodis
,
M.
Chachisvilis
,
B. R.
Cherry
,
M.
Fournier
,
A.
Moore
,
T. A.
Moore
, and
D.
Gust
, “
Multiporphyrin arrays with π–π interchromophore interactions
,”
J. Am. Chem. Soc.
137
,
245
258
(
2015
).
26.
L.
Moretti
,
B.
Kudisch
,
Y.
Terazono
,
A. L.
Moore
,
T. A.
Moore
,
D.
Gust
,
G.
Cerullo
,
G. D.
Scholes
, and
M.
Maiuri
, “
Ultrafast dynamics of nonrigid zinc-porphyrin arrays mimicking the photosynthetic ‘special pair
,’”
J. Phys. Chem. Lett.
11
,
3443
3450
(
2020
).
27.
Y.
Terazono
,
E. J.
North
,
A. L.
Moore
,
T. A.
Moore
, and
D.
Gust
, “
Base-catalyzed direct conversion of dipyrromethanes to 1,9-dicarbinols: A [2 + 2] approach for porphyrins
,”
Org. Lett.
14
,
1776
1779
(
2012
).
28.
S.
Cho
,
W.-S.
Li
,
M.-C.
Yoon
,
T. K.
Ahn
,
D.-L.
Jiang
,
J.
Kim
,
T.
Aida
, and
D.
Kim
, “
Relationship between incoherent excitation energy migration processes and molecular structures in zinc(II) porphyrin dendrimers
,”
Chem. - Eur. J.
12
,
7576
7584
(
2006
).
29.
G.
Kodis
,
Y.
Terazono
,
P. A.
Liddell
,
J.
Andréasson
,
V.
Garg
,
M.
Hambourger
,
T. A.
Moore
,
A. L.
Moore
, and
D.
Gust
, “
Energy and photoinduced electron transfer in a wheel-shaped artificial photosynthetic antenna-reaction center complex
,”
J. Am. Chem. Soc.
128
,
1818
1827
(
2006
).
30.
M.
Kretov
,
I.
Iskandarova
,
B.
Potapkin
,
A.
Scherbinin
,
A.
Srivastava
, and
N.
Stepanov
, “
Simulation of structured 4T16A1 emission bands of Mn2+ impurity in Zn2SiO4: A first-principle methodology
,”
J. Lumin.
132
,
2143
2150
(
2012
).
31.
M.
Kretov
,
A.
Scherbinin
, and
N.
Stepanov
, “
Simulating the structureless emission bands of Mn2+ ions in ZnCo3 and CaCo3 matrices by means of quantum chemistry
,”
Russ. J. Phys. Chem. A
87
,
245
251
(
2013
).
32.
P. V.
Yurenev
,
M. K.
Kretov
,
A. V.
Scherbinin
, and
N. F.
Stepanov
, “
Environmental broadening of the CTTS bands: The hexaammineruthenium(II) complex in aqueous solution
,”
J. Phys. Chem. A
114
,
12804
12812
(
2010
).
33.
Z.
Shuai
,
H.
Geng
,
W.
Xu
,
Y.
Liao
, and
J.-M.
André
, “
From charge transport parameters to charge mobility in organic semiconductors through multiscale simulation
,”
Chem. Soc. Rev.
43
,
2662
2679
(
2014
).
34.
P. S.
Rukin
,
A. Y.
Freidzon
,
A. V.
Scherbinin
,
V. A.
Sazhnikov
,
A. A.
Bagaturyants
, and
M. V.
Alfimov
, “
Vibronic bandshape of the absorption spectra of dibenzoylmethanatoboron difluoride derivatives: Analysis based on ab initio calculations
,”
Phys. Chem. Chem. Phys.
17
,
16997
17006
(
2015
).
35.
G. J.
Small
, “
Herzberg–Teller vibronic coupling and the Duschinsky effect
,”
J. Chem. Phys.
54
,
3300
3306
(
1971
).
36.
M.
Lax
, “
The Franck-Condon principle and its application to crystals
,”
J. Chem. Phys.
20
,
1752
1760
(
1952
).
37.
F.
Neese
,
T.
Petrenko
,
D.
Ganyushin
, and
G.
Olbrich
, “
Advanced aspects of ab initio theoretical optical spectroscopy of transition metal complexes: Multiplets, spin-orbit coupling and resonance Raman intensities
,”
Coord. Chem. Rev.
251
,
288
327
(
2007
).
38.
M. D.
Frank-Kamenetskiĭ
and
A. V.
Lukashin
, “
Electron-vibrational interactions in polyatomic molecules
,”
Sov. Phys. - Usp.
18
,
391
409
(
1975
).
39.
D.-Y.
Wu
,
M.
Hayashi
,
Y.-J.
Shiu
,
K.-K.
Liang
,
C.-H.
Chang
, and
S.-H.
Lin
, “
Theoretical calculations on vibrational frequencies and absorption spectra of S1 and S2 states of pyridine
,”
J. Chin. Chem. Soc.
50
,
735
744
(
2003
).
40.
C. A.
Rozzi
,
F.
Troiani
, and
I.
Tavernelli
, “
Quantum modeling of ultrafast photoinduced charge separation
,”
J. Phys.: Condens. Matter
30
,
013002
(
2017
).
41.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
D.
Williams-Young
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
,
Gaussian ∼ 16 Revision C.01
,
Gaussian, Inc.
,
Wallingford, CT
,
2016
.
42.
T.
Yanai
,
D.
Tew
, and
N.
Handy
, “
A new hybrid exchange–correlation functional using the Coulomb-attenuating method (CAM-B3LYP)
,”
Chem. Phys. Lett.
393
,
51
57
(
2004
).
43.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
, “
Quantum mechanical continuum solvation models
,”
Chem. Rev.
105
,
2999
3094
(
2005
).
44.
R. L.
Martin
, “
Natural transition orbitals
,”
J. Chem. Phys.
118
,
4775
4777
(
2003
).
45.
T.
Lu
and
F.
Chen
, “
Multiwfn: A multifunctional wavefunction analyzer
,”
J. Comput. Chem.
33
,
580
592
(
2012
).
46.
P. J.
Spellane
,
M.
Gouterman
,
A.
Antipas
,
S.
Kim
, and
Y. C.
Liu
, “
Porphyrins. 40. Electronic spectra and four-orbital energies of free-base, zinc, copper, and palladium tetrakis(perfluorophenyl)porphyrins
,”
Inorg. Chem.
19
,
386
391
(
1980
).
47.
F.
Santoro
,
A.
Lami
,
R.
Improta
,
J.
Bloino
, and
V.
Barone
, “
Effective method for the computation of optical spectra of large molecules at finite temperature including the Duschinsky and Herzberg–Teller effect: The Qx band of porphyrin as a case study
,”
J. Chem. Phys.
128
,
224311
(
2008
).
48.
B.
Minaev
,
Y.-H.
Wang
,
C.-K.
Wang
,
Y.
Luo
, and
H.
Ågren
, “
Density functional theory study of vibronic structure of the first absorption Qx band in free-base porphin
,”
Spectrochim. Acta, Part A
65
,
308
323
(
2006
).
Published open access through an agreement with Istituto Nanoscienze Consiglio Nazionale delle Ricerche

Supplementary Material