The electronic structure of eight zinc-centered porphyrin macrocyclic molecules are investigated using density functional theory for ground-state properties, time-dependent density functional theory (TDDFT) for excited states, and Franck-Condon (FC) analysis for further characterization of the UV-vis spectrum. Symmetry breaking was utilized to find the lowest energy of the excited states for many states in the spectra. To confirm the theoretical modeling, the spectroscopic result from zinc phthalocyanine (ZnPc) is used to compare to the TDDFT and FC result. After confirmation of the modeling, five more planar molecules are investigated: zinc tetrabenzoporphyrin (ZnTBP), zinc tetrabenzomonoazaporphyrin (ZnTBMAP), zinc tetrabenzo*cis*diazaporphyrin (ZnTB*cis*DAP), zinc tetrabenzo*trans*diazaporphyrin (ZnTB*trans*DAP), and zinc tetrabenzotriazaporphyrin (ZnTBTrAP). The two latter molecules are then compared to their phenylated sister molecules: zinc monophenyltetrabenzotriazaporphyrin (ZnMPTBTrAP) and zinc diphenyltetrabenzo*trans*diazaporphyrin (ZnDPTB*trans*DAP). The spectroscopic results from the synthesis of ZnMPTBTrAP and ZnDPTB*trans*DAP are then compared to their theoretical models and non-phenylated pairs. While the Franck-Condon results were not as illuminating for every B-band, the Q-band results were successful in all eight molecules, with a considerable amount of spectral analysis in the range of interest between 300 and 750 nm. The *π*-*π*^{∗} transitions are evident in the results for all of the Q bands, while satellite vibrations are also visible in the spectra. In particular, this investigation finds that, while ZnPc has a *D*_{4h} symmetry at ground state, a *C*_{4v} symmetry is predicted in the excited-state Q band region. The theoretical results for ZnPc found an excitation energy at the Q-band 0-0 transition of 1.88 eV in vacuum, which is in remarkable agreement with published gas-phase spectroscopy, as well as our own results of ZnPc in solution with Tetrahydrofuran that are provided in this paper.

## I. INTRODUCTION

Properties of metallophthalocyanines (MPcs) have been investigated for nearly a century, beginning with a series of articles published by Linstead and Coworkers in the early 1930s.^{1–3} Since then, MPcs have been studied for use in chemical sensors, photodynamic therapy, optics, bioelectronics, molecular optoelectronic gates, and molecular neural networks.^{4–7} Of particular importance is their application in organic semiconductors, discovered in 1948 by Eley.^{8} Metal-centered phthalocyanines were first introduced in a laminate in 1958, but not until 1986 was the first report of a heterojunction organic photovoltaic (OPV) cell published.^{9,10} Since Tang demonstrated the first efficient OPV at 1% power conversion efficiencies (PCEs), phthalocyanines have been a mainstay in the organic electronics field.^{10} Single heterojunction devices employing MPcs are now surpassing 6% PCE with a modified zinc phthalocyanine material.^{11–13} The targeted synthesis of selected materials based on computational models would greatly serve to expedite the identification of optimal materials and accelerate progress in organic electronics as well as the wide reaching fields of chemistry and physics. The present research thus utilizes a theoretical approach to investigate zinc phthalocyanine (ZnPc), zinc tetrabenzoporphyrin (ZnTBP), and azaporphyrin analogs. Fig. 1 shows a diagram of ZnPc where X = N, ZnTBP where X = C−H, and the porphyrin analogs will have varying N, C–H, or Phenyl bonds to X in the *meso* position.

UV spectra of metallophthalocyanines and metallotetrabenzoporphyrins were first given spectral assignments by Edwards, Gouterman, and co-workers.^{14,15} These earlier band assignments were based on the Gouterman four-orbital model.^{14–16} There have been ZnPc studies proposing a symmetry forbidden *n*-*π*^{∗} transition blue-shifted from the highest energy states of the Q band region.^{17–20} Even though the additional assignments in the first four bands have been debated in literature, Gouterman and Edwards made a key observation in 1970. They noted that the Q band, in particular, would be subject to Franck-Condon (FC)^{21–23} vibrational displacement, and point to an earlier publication by Gouterman and Fulton, in which the higher energy bands were believed to take on much more peculiar shapes, especially at the investigated temperature of 0 K.^{24} Some believe the four orbital model is simplistic, but the framework is generally still reliable for the major transitions, and is used in this present work.

Most functionals used for time-dependent density functional theory (TDDFT) have not been extremely accurate, leading to low energy excitations that are usually blue shifted. For higher excitation energies, more states must be resolved, which may be blue or red shifted. The first TDDFT study of ZnPc by Ricciardi *et al.* used an adiabatic local-density approximation (ALDA) in combination with a statistical average of different orbital potentials (SAOP) for the exchange-correlation potential to solve for the time-dependent density.^{25} The Q-band result for the 0-0 transition was close to experimental values: 1.96 eV.^{25} The significant contribution of Ricciardi *et al.* was to separate and change the exchange-correlation energy. It has subsequently become much more clear that the accuracy of excited-state energies is more dependent on the percentage of Hartree-Fock (HF) exchange or exchange energy.^{18,25,26} The current research explores these unique molecules as seen is Fig. 2 by looking at range-separated hybrid functionals for excited states.

Theoretical research has historically studied ZnPc and ZnTBP more often along with porphyrin analogs that have a *D*_{4h} symmetry.^{18,25,27} Recently, these molecules along with a copper center and asymmetric analogs have been synthesized and investigated with magnetic circular dichroism (MCD) theory.^{28} In fact, there is a large amount of investigation of MPcs over the last couple of decades by Mack, Stillman, and Kobayashi.^{19,20,28–30} However, there has been minimal research presented on these molecules using Franck-Condon analysis until a recent study by Gou *et al.* presented work on ZnPc and ZnTBP.^{27} Two approaches seem to have dominated the theoretical work. One uses a combination of density functional theory (DFT), TDDFT, and oscillator strengths (intensities) of excitation energies. The other uses MCD theory, as well as some DFT. This current research employs techniques and theoretical framework from the former using DFT and vertical excitation energies found with TDDFT. In addition, this research uses further optimization of each excited state leading to better approximated intensities and incorporating the FC analysis. The FC analysis can become computationally expensive due to the necessary TDDFT calculations for frequency, while there are increasing difficulties in resolving excited-state energies with TDDFT as they approach the B band and higher.

The specific molecules that have been studied here are zinc-centered and consist of zinc phthalocyanine (ZnPc), zinc tetrabenzoporphyrin (ZnTBP), zinc tetrabenzomonoazaporphyrin (ZnTBMAP), zinc tetrabenzo*cis*diazaporphyrin (ZnTB*cis*DAP), zinc tetrabenzo*trans*diazaporphyrin (ZnTB*trans*DAP), zinc diphenyltetrabenzo*trans*diazaporphyrin (ZnDPTB*trans*DAP), zinc tetrabenzotriazaporphyrin (ZnTBTrAP), and monophenyltetrabenzotriazaporphyrin (ZnMPTBTrAP), as can be seen in Fig. 2. The zinc-centered molecules are the focus due to their superior performance as a leading MPc for absorption in OPVs amongst other possible metal centers.^{31} Furthermore, the research into phenyl ring additions is due to the previously correlated research where the increase in steric bulk leads to better PCE.^{13,32} This current research investigates all 8 of the molecules in Fig. 2 for their excited states and absorption in the UV-vis region using TDDFT and FC analysis. This class of molecules required symmetry breaking in several of the excited states to find accurate geometries in order to incorporate the FC analysis, as will be discussed in Sec. II.

## II. THEORETICAL AND COMPUTATIONAL METHODS

All of the computational work in this paper used the Gaussian 09 software and GaussView 5.^{33} All the calculations were done on Trestles, a dedicated Extreme Science and Engineering Discovery Environment (XSEDE) cluster designed by Appro and San Diego Supercomputer Center (SDSC). The ground-state calculations in this paper were all done based on the Kohn-Sham DFT.^{34} The Becke, three-parameter, Lee-Yang-Parr (B3LYP) hybrid functional was used for all the Molecular Orbital’s (MO’s) energies to screen the HOMO LUMO levels of each molecule.^{35,36} B3LYP was benchmarked against several other methods, including the M11^{37} functional that is used for the rest of the calculations. The results show that M11 overestimated the MO levels as compared to most other functionals; this can be found in the supplementary material.^{38} This is why the B3LYP method was only used for the MO levels. The B3LYP functional had already been reliably proven by several researchers to give excellent agreement with experiment for ground-state properties of ZnPc.^{18,26} However, all the ground-state molecular geometry optimization and frequency calculations were done with the M11 functional. The ground-state structure is benchmarked with several functionals in the supplementary material^{38} and the results shows very little difference in the bond calculations.

The excited-state calculations were done based on TDDFT, an extension of the Kohn-Sham equations.^{39,40} The excited-state research uses a range-separated hybrid meta-generalized gradient approximations (GGA) method called M11 from the Truhlar group, in which a HF exchange can range from 42.8% to 100% (see Ref. 37 for the general formula for this exchange model). The M11 functional was chosen based on benchmarking for excited-state methods using ZnPc as the baseline; see Table I. Only the initial time-dependent energy is given in Table I, the optimization of each excited state is a separate matter. While many of these functionals overestimate the energies as compared to experiment, a few of them actually underestimate the energy in the B band. These options did not seem reliable for the purpose since the energy will become lower as the excited state is optimized. In the initial results, M11 was chosen for having the lowest energy for the 0-0 transition in the Q band as compared to experiment. The only other range separated functional considered was Coulomb-attenuating method (CAM)-B3LYP.^{41} While CAM-B3LYP looked good over all for both bands, M11 was chosen based on the accuracy for the Q-band 0-0 transition.

Method . | 6-31g(d) (eV) . | 6-31+g(d) (eV) . | Expt.^{a} (eV)
. |
---|---|---|---|

B3LYP | 2.09 | 2.05 | 1.88 |

3.38 | 3.35 | 3.80 | |

PBE0 | 2.13 | 2.10 | 1.88 |

3.55 | 3.52 | 3.80 | |

CAM-B3LYP | 2.07 | 2.03 | 1.88 |

4.02 | 3.97 | 3.80 | |

M062X | 2.14 | 2.10 | 1.88 |

3.97 | 3.97 | 3.80 | |

M06 | 2.04 | 2.00 | 1.88 |

3.51 | 3.46 | 3.80 | |

M11 | 2.01 | 1.96 | 1.88 |

4.27 | 4.22 | 3.80 |

Method . | 6-31g(d) (eV) . | 6-31+g(d) (eV) . | Expt.^{a} (eV)
. |
---|---|---|---|

B3LYP | 2.09 | 2.05 | 1.88 |

3.38 | 3.35 | 3.80 | |

PBE0 | 2.13 | 2.10 | 1.88 |

3.55 | 3.52 | 3.80 | |

CAM-B3LYP | 2.07 | 2.03 | 1.88 |

4.02 | 3.97 | 3.80 | |

M062X | 2.14 | 2.10 | 1.88 |

3.97 | 3.97 | 3.80 | |

M06 | 2.04 | 2.00 | 1.88 |

3.51 | 3.46 | 3.80 | |

M11 | 2.01 | 1.96 | 1.88 |

4.27 | 4.22 | 3.80 |

^{a}

Observed data from Ref. 14.

While several basis sets were considered, 6-31g(d) was used throughout.^{42–44} Both 6-31g(d) and the 6-31+g(d) basis sets were seriously considered from the initial benchmarking. However, the extra diffuse function in the 6-31+g(d) basis set turned out to be problematic for the convergence of the TDDFT calculations, and only yielded negligible bond energy differences at the ground state. There was a notable difference in the excited-state results, as can be seen in Table I, but the small difference was not worth the computational cost and efficiency (this can be seen in Tables II and III in the supplementary material^{38}).

The TDDFT energy results were then used to optimize the TDDFT geometry of each excited state separately in the UV-vis range of interest. Once each excited state was correctly converged on, a TDDFT calculation for frequency was then carried out. These calculations are not trivial, as they require computing second derivatives of the energy obtained by numerical differentiation of the analytic first derivatives. The results of these calculations are in the form of the diagonalization of the force constants matrix. This produces either positive or negative eigenvalues, where the negative eigenvalues correspond to imaginary frequencies. This outcome actually gives us much needed information. Each of the imaginary frequency cases indicated a saddle point of order 1, transition state, or a possible higher order on the potential energy surface (PES). The molecule has to be re-optimized by breaking the symmetry along the particular frequency mode with the largest negative eigenvalue. This displacement allows the molecule to take on a non-planar form for the optimization, which was required for several of the excited states in this research.

Once the TDDFT frequency converged correctly, the FC analysis could produce the spectrum by using the vibrational data from the ground-state and excited-state frequencies. The exception occurs in cases with no overlap or when the progression of the overlap integrals are too low between the initial (ground) and the final (excited) states. The FC principle evolved out of the theory by Franck to estimate vibrational electronic intensities for transitions in diatomic molecules.^{23,45,46} Condon soon developed an expansion of this idea with a way to calculate probability amplitudes for electronic vibrational transitions, given by the Condon overlap integral in Eq. (1). The original formalism only considered a harmonic oscillator solution. The *R* represents the nuclear interdistance, *n* is the vibrational quantum number of the excited state, and *n*′ is the quantum number of the ground state,^{21,47}

This method further assumes the ability to separate the electronic, rotational, and vibrational parts of the wave functions. In this solution, the transition probabilities are proportional to the square of the vibrational overlap integrals *C* between the final and initial state.^{48} The overlap integrals are considered the Franck-Condon factors (FCFs). There have been many historical explanations of Franck-Condon factors.^{48–50} Performing such first-principles calculations is challenging because each excited state needs to be optimized separately, which can produce mode mixing first described by Duschinsky.^{51} In some cases, after the calculation of the TDDFT frequency, the solution produces imaginary frequency solutions, which indicates that the optimization for that excited state did not find the true minimum on the PES. One must then employ symmetry breaking to find the true minimum, and in turn the correct minimum on the PES, in order for the FC analysis to calculate the correct transition.

Once the framework for the FC integrals is set up, there could be concern over the ability of the Cartesian coordinate system to accurately represent distortions of the bonds correctly due to coupling in the Duschinsky matrix. This concern is rooted in whether the weighted Cartesian displacement vectors can accurately express the shift in coordinates for the computed PES. This has been debated in recent years between Cartesian and internal coordinate systems.^{52,53} While it has been shown that internal coordinate system yield better results for the vibronic equations due to the fact that there are less coordinates displaced versus the Cartesian system, some of the same researchers have produced significant work and comparisons in Cartesian coordinates.^{54,55} In addition to comparisons between the adiabatic and vertical Hessian.^{56} It should be noted that an internal coordinate system is many times preferred, but not as widely available. This work presented here has been carried out in Cartesian coordinates. The results are satisfactory as compared to the experimental data that were available to this research. As far as any doubts about the Cartesian coordinate system, ZnPc is generally a relatively rigid molecule. Further evidence of this is discussed in Sec. V A, where there are quantitative results given for the small energy difference between the optimized geometry of the excited state where an imaginary frequency occurs (the saddle point) and the true minimum geometry found with the broken symmetry of ZnPc.

There are several approaches in literature that give in-depth discussions on these calculations.^{50,51,57–60} A few key equations from the harmonic solution are given here, but the entire solution is calculable, and can be found in Santoro *et al.*^{57} The basic idea is that the molar absorption is produced by the absorption equation, which gives the transition between two electronic states |Ψ_{w′}〉 and |Ψ_{w}〉,^{61}

This profile in Eq. (2) is often referred to as the *stick absorption* and will be seen in many TDDFT results with a single vertical line on the spectrum. This is often fitted with a gaussian broadening distribution function. The molecular initial and final wave functions are Ψ_{w′} and Ψ_{w}, respectively, the transition dipole moment is represented by *μ*, and the intensity of the absorption is governed by the 〈Ψ_{w′}|*μ*|Ψ_{w}〉 integral. The electronic transition dipole moment can be expressed in a power series of the normal coordinates **Q**,^{57}

Here, the *N*-dimensional vibrational states of |w′〉 and |w〉 are products of |w_{k}〉 for each mode *k*, |w〉 ≡ |**w**〉 = |w_{1}〉 ⊗ |w_{2}〉⋯ ⊗ |w_{N}〉.^{57} We can obtain the overlap integrals 〈**w**′|**w**〉. Where **Q** represents the normal coordinates for the initial state and **Q**′ the normal coordinates for the final state. The transformation matrix is represented by **J** and the equilibrium displacement vector **K**.^{57,58} This treatment is expanded, using the Duschinsky transformation to obtain a common coordinate set,^{49}

While one significant improvement to the calculations came when Doktorov used Glauber’s coherent state to solve for the FCFs, another solution came in the form of generating functions from Sharp and Rosenstock’s work.^{48,60,62} This allows the Gaussian program to compute the spectrum band-shape as a sum of individual transitions as shown in Eq. (2) by means of recursion formulas. This application is used to find the general overlap *I*_{w′,w} = 〈**w**′|**w**〉.^{57} Each transition is broadened by means of the default distribution Gaussian function in the Gaussian 09^{33} program using a half width half max (HWHM) of 135.00 cm^{−1}.

## III. EXPERIMENTAL SECTION

Zinc phthalocyanine was purchased from Sigma Aldrich and ZnDPTBtransDAP and ZnMPTBTrAP were synthesized according to previous reported procedure.^{63} All materials were purified via sublimation in a four zone thermal gradient furnace prior to use. The UV-visible spectra were recorded on a Hewlett-Packard 4853 diode array spectrometer in a dilute solution with tetrahydrofuran (THF) as the solvent.

## IV. MOLECULAR ORBITAL LEVELS

The molecular orbital calculations were all done with the B3LYP functional as well as the 6-31g(d) basis set. The first set of planar porphyrin analog results is in Fig. 3, the top half of the figure has the calculated 6 states closest to the HOMO LUMO gap. The extra HOMO states for ZnPc were calculated as well to complete the *a*_{2u} four orbital assignment. The major orbital assignments are labeled in a dotted line across Fig. 3. The bottom half of the figure shows each of these molecular orbitals corresponding to their respective orbital energies in the top half of Fig. 3. Notice that the orbitals with degeneracies have exact overlaps in the LUMO energies, such as ZnPc and ZnTBP. The LUMO orbitals that are degenerate or nearly degenerate have identical orbitals with an x orientation or y orientation, as expected. If you look at a molecule such as ZnTBTrAP, there is clearly a separation in the energy of the LUMO and the LUMO + 1 state, this can be attributed to the asymmetry and position of the nitrogens, which has already been reported.^{64} The numerous orbital configurations that were calculated are not all included here, since the focus of this research encompasses 8 molecules in total. Another energy level to point out is the energy level of the LUMO and LUMO + 1 for ZnTB*cis*DAP, these two energies are so close (right on top of each other) in energy, they are nearly degenerate. The four orbital model sits in the HOMO-1 position for all of the azaporphyrins investigated, except for ZnPc, where the four orbitals are found at the HOMO-5 position.

From Fig. 3, it is apparent that there are two molecules with lower HOMO LUMO gaps out of the rest, as compared to ZnPc. These include ZnTB*trans*DAP and ZnTBTrAP. The results for these two molecules find the four orbital model, as described by Gouterman and published by others for transition state contributions, in the HOMO-1 position.^{18,65} This orbital is considered the minor contribution to the LUMO level, the fact that these two molecules have this orbital sitting closer to the HOMO level as compared to ZnPc may be another reason to take a closer look at these two for their transitions as compared to ZnPc. While ZnPc is already well investigated as a good donor molecule, the present research looks at the electronic structure and the theoretical absorption of all of these porphyrin analogs to see if there are other promising donor molecules.

The addition of the phenyl rings was done to find a molecule that has a combination of properties such as a low HOMO LUMO gap, while having the advantage of lending some steric bulk to the molecule in the prediction that it will help prevent recombination of excitons created after absorption into an OPV.^{13,32} ZnTBTrAP and ZnTB*trans*DAP were clearly two molecules of interest that came out of the molecular orbital levels as seen in Fig. 3 where the HOMO LUMO gaps were 2.20 eV and 2.19 eV, respectively.

Ground-state optimizations with M11 and 6-31g(d) were done for all the azaporphyrin analogs with phenyl additions, only two stayed relatively planar keeping the phenyl addition orthogonal, whereas the other four had considerable buckling and saddle shaped ground states. These were coincidentally the same two molecules that had the lowest HOMO LUMO gap. The two molecules, zinc monophenyltetrabenzotriazaporphyrin (ZnMPTBTrAP) and zinc diphenyltetrabenzo*trans*diazaporphyrin (ZnDPTB*trans*DAP) can be seen in Fig. 2. Their comparison molecular orbitals and energy levels can be seen in Fig. 4. The additional phenyl groups decreased the HOMO LUMO gap compared to their original analogs. Both ZnDPTB*trans*DAP and ZnMPTBTrAP have results of lower HOMO LUMO gaps, 2.16 eV and 2.18 eV, respectively, compared to the ZnPc gap of 2.19 eV. This comparison in Fig. 4 now makes these two molecules of high interest to investigate further. This was a desired outcome, since there might be an advantage of intermolecular interactions with the planarity in the inner macrocycle as well as the added steric bulk.

## V. RESULTS AND DISCUSSION

### A. ZnPc: Results and benchmarking

The following benchmarking was done with ZnPc, since the excited states are well-characterized. Results appearing in Table I show better agreement was achieved for the 0-0 transition of the Q-band edge with the 6-31+g(d) basis set rather than the 6-31g(d) basis set. However, the extra diffuse function resulted in much longer calculation times, and in a couple of cases, the frequencies did not converge with the TDDFT level of theory. This was due to the diffuse function creating additional linearly dependent equations on top of the numerical second derivatives that are calculated in the TDDFT frequency calculations. Benchmarking was done again to find a method that would work with the 6-31g(d) basis set. The M11 method was decided on because it came closest to the Q-band 0-0 transition, even though it was a bit worse for the B-band region than other possibilities.

As discussed in Sec. II, the two range separated methods, M11 and CAM-B3LYP, were two of the methods that were of high interest. CAM-B3LYP did show promising results for a more balanced UV-vis prediction, but both methods over estimated the 0-0 transition of the B-band edge. The M11 method was chosen for the accuracy in the TD approximation to the 0-0 transition of the Q-band edge. The 0-0 transition energies were not corrected for zero point vibrational energy (ZPVE) differences. In addition, it is fair to assume that the difference in ZPVEs between the two states is likely to be small. However, it is important to note that in the final Franck-Condon analysis, the ZPVEs are automatically taken into account in the 0-0 transition as well as all of the transitions.

The next step was taking each TD energy and running a TD optimization on each state to confirm the correct geometry in the excited state. These calculations optimized the TD geometry, lowered the energy on the new PES, and in several degenerate cases it would optimize for one out of the two degenerate states, since only one can be converged on. This results in the vertical excitation solutions or “stick absorption” seen in Fig. 5 and the second column of Table II.

State . | TDDFT energy . | f
. | Assignment . | Symmetry . | Experiment (eV) . | ||
---|---|---|---|---|---|---|---|

ZnPc | |||||||

1A_{1g} | Ground state | D_{4h} | |||||

1B_{2} | 1.88 eV (659 nm) | 0.5318 | Q (π-π∗) | C_{4v} → C_{2v} | 1.88^{a} | 1.89^{b} | 1.86^{c} |

1B_{1} | 2.10 eV (590 nm) | 0.4449 | Q (π-π∗) | C_{4v} → C_{2v} | 2.07^{a} | 2.07^{c} | |

1A | 3.29 eV (377 nm) | 0.0960 | B (π-π∗) | C_{2} → C_{2} | 3.80^{a} | 3.71^{b} | 3.63^{c} |

2A | 3.82 eV (324 nm) | 0.0922 | B (π-π∗) | C_{2} → C_{2} | 3.74^{b} |

State . | TDDFT energy . | f
. | Assignment . | Symmetry . | Experiment (eV) . | ||
---|---|---|---|---|---|---|---|

ZnPc | |||||||

1A_{1g} | Ground state | D_{4h} | |||||

1B_{2} | 1.88 eV (659 nm) | 0.5318 | Q (π-π∗) | C_{4v} → C_{2v} | 1.88^{a} | 1.89^{b} | 1.86^{c} |

1B_{1} | 2.10 eV (590 nm) | 0.4449 | Q (π-π∗) | C_{4v} → C_{2v} | 2.07^{a} | 2.07^{c} | |

1A | 3.29 eV (377 nm) | 0.0960 | B (π-π∗) | C_{2} → C_{2} | 3.80^{a} | 3.71^{b} | 3.63^{c} |

2A | 3.82 eV (324 nm) | 0.0922 | B (π-π∗) | C_{2} → C_{2} | 3.74^{b} |

In the case of ZnPc, the symmetry had to be broken in order to achieve the lowest energy, as described in Sec. II. The symmetry for the lowest energy of each excited state is listed in column 5 of Table II. The result of a *C*_{4v} symmetry was found for the Q-band excited state, which was predicted by Scheidt and Dow based on the Zn–N bond distances, where the Zn center could stick out of the plane.^{66} This conclusion was reported by two additional groups in the 1990s based on the idea that this could be observed in a deoxygenated state, in addition to possible Jahn-Teller splitting.^{67,68} In fact, the displacement carried out along the normal coordinates associated to the imaginary frequency moved the zinc center out of the plane by 0.24 Å, which happens to be in close agreement to experiment findings by Mihill *et al.*^{67} The present theoretical results show a *D*_{4h} symmetry is preserved in the ground state, but the symmetry breaks in the first and second degenerate excited states. The Q-band excited state could be attributed to a small split, or nearly degenerate excited state, but these results show the B band to be a much larger split. In addition, the result of the 0-0 transition in the Q band was 1.88 eV, exactly the same energy as observed in gas phase, as can be seen in Table II. Since these calculations were done in vacuum, gas phase would be the closest to compare with. The experimental data in Fig. 5 were taken in THF solvent. There would be an expected blue-shift of gas phase data as compared to solution phase data such as the spectrum given in Fig. 5.

The TDDFT energy “stick absorption” results as discussed in Sec. II are presented in Table II and shown in Fig. 5. The Franck-Condon analysis results from the TDDFT frequency calculations are also shown in Fig. 5, as well as the present spectrum in THF (black solid line). The FC results for the Q band produced a spectrum (red dotted line) that is remarkably similar to the experimental spectrum. The FC spectrum is shifted from the experimental spectrum, one would expect a slight blue shift from spectral analysis taken in solution. In addition to the fact that the DFT and TDDFT frequency calculations were done at 0 K and the spectroscopy was done at 298 K. There are small peaks that can be seen at the higher energy end of the Q band where there could be vibrational transition satellites. One might conclude this as evidence of a *n*-*π*^{∗} symmetry forbidden transition, but vibronically allowed coupled state. This result can only be observed qualitatively as an outcome, but cannot be quantified here, since the result is from an adiabatic transition. Despite the limitations of using Cartesian coordinates and an adiabatic transition, it is quite remarkable that there are peaks that resemble symmetry forbidden peaks, in addition to the *π*-*π*^{∗} absorption energies found so close to experiment as seen in Table II.

The result of the first degenerate state breaks into two states, this has been speculated by others both in theory and in experiment.^{20,25,68,69} This research is reported as having surprising results worth further study for absorption using FC analysis. Even though the beginning state would produce an *E _{u}* transition state as is widely reported, these results give two states that are

*B*

_{2}and

*B*

_{1}. With a starting symmetry of

*D*

_{4h}→

*C*

_{4v}, hence it is not surprising that the molecule converged with no negative eigenvalues at a

*C*

_{2v}, since there is no doublet representation for

*C*

_{4v}, this is a plausible outcome.

There is a clarification to the results of the Q band to be addressed here. One might conclude that due to the removal of symmetric coordinates a solution of a double well could be produced in the excited state PES. This is a serious concern, in any case, when one starts at a high degree of symmetry such as the *D*_{4h} in ZnPc, the removal of symmetry to a subgroup by one degree indicates a double well. In the case that the symmetry could be broken to a 2nd degree subgroup, this can indicate multiple wells in the PES. The main concern lies between the minimum and just above the barrier (or saddle) in the PES. There are potential quantum mechanical effects that could be due to tunneling near the edge of the barrier and distortion of the wave function if the next excited state lied just above the barrier. There is an easy way to check whether the removal of symmetry is problematic. The original energy before the symmetry is broken is where the saddle or barrier would reside. In the case of the Q band in ZnPc, the energy difference between the minimum and the barrier is 0.0023 eV. If one considers doubling this energy as a safe buffer from the barrier, then the next excited state has to come 0.0046 eV above the minimum. Since the next excited state in Table II is at 2.10 eV, there is a large enough energy buffer not to be effected by the barrier.

In addition, the global minimum in ZnPc has only a small energy difference (0.002 eV) compared to the saddle point on the PES, so the barrier is effectively a perturbation. So, the use of the harmonic approximation for the anharmonic double-well, although not exact, is reasonable. Also, one can speculate that the coupling to the normal mode is stronger than to the mode at the saddle point.

Although the ZnPc transition for the Q band had excellent agreement with experiment, there are excitations where time-dependent theory and the spectral analysis do not agree. The B band has a TDDFT energy that is underestimated in energy as compared to the experimental data and the original TD benchmark in Table I. In the benchmarking, the initial energy was predicted at nearly 1 eV higher than the end result after the TD optimization of the 2nd degenerate excited state of ZnPc, and the result was found with a *C*_{2} symmetry for the lowest energy. This appeared to be a bad result at first, but one has to consider the possibility that the two states found in the B band may not have been degenerate, and the two energies do fall under the broad absorption of the B Band. In fact, these two states come very close to the energies that Nyokong *et al.* found in the B band of ZnPc in cyanide at 3.21 eV (386 nm) and 3.75 eV (331 nm).^{70} The present calculation finds two states at 3.29 eV (377 nm) and 3.82 eV (324 nm), since the symmetry was broken to *C*_{2}, both of these energies are found as A states, one would expect a B state as a representation of the original doublet state.

The Frank-Condon analysis did not produce much of a spectral analysis in the B band, this was because the FC analysis was only able to progress to 70% of the overlap integrals. This could be due to the large energy shift and a very diffuse matrix. There is also evidence in the B band that there could be an issue with multiple wells in the excited state PES. Here, the difference between the barrier (or saddle) and the minimum has a much larger energy difference at 0.458 eV. This energy difference cannot be attributed to a perturbation and the FC resolution of the B band cannot be considered reliable. For this reason and the minimal resolution for the B band in ZnPc, the rest of the results will only focus on the FC analysis for the Q bands while including TDDFT optimized energy results throughout. However, the Q-band results were encouraging enough to move forward on the rest of the azaporphyrin molecules, where there is minimal experimental data for comparison.

### B. ZnTBP and azaporphyrin analog results

Before discussion of this set of molecules, there is an important theoretical point to be clear about. Since the ZnTBP molecules has a *D*_{4h} symmetry and is doubly degenerate in the excited state, there was only one state to converge on for each of the Q band and the B band. The azaporphyrin analogs, on the other hand, had less symmetry to begin with, therefore, no less than 4 excited states needed to be calculated (two for each the Q and B band region). Because of the unique absorption of these zinc phthalocyanine and porphyrin analogs, the Q and B bands all had a degree of energy separation, meaning there is not much absorption or “noise” between bands, this required a methodology or protocol for solving the FC analysis. Each band had to be resolved by using the lowest energy excitation of the band edge. The FC results would envelope the whole band and resolve any other excited states within that band.

There is evidence in this section that ZnTBP proves to be a good benchmark for the 0-0 transition of the Q-band edge using the M11 functional as well. In addition to ZnTBP, the zinc-centered azaporphyrin series that was investigated was a set of planar molecules with altered macrocyclic structures in the *meso* position ranging in degree of symmetry and nitrogens: ZnTBMAP, ZnTB*cis*DAP, ZnTB*trans*DAP, and ZnTBTrAP as compared to ZnPc. These can be seen as the non-phenylated 6 molecules in Fig. 2. These azaporphyrins were first given assignments according to Gouterman’s four orbital model.^{65}

The excited-state results are then investigated and the absorption of all the azaporphyrin analog molecules are looked at with the same methodology described in Sec. V A for ZnPc. The rest of the porphyrin analogs were optimized using M11 and 6-31g(d) from the ground state, calculated for the time-dependent states using TDDFT, converged on for each TD state, and if there was symmetry breaking required, it is all described in Table III. The TD optimized states or “stick absorption” are shown in Fig. 6. The theoretical Franck-Condon absorption results are also overlaid in Fig. 6 as well. There was Franck-Condon convergence for every Q band in each of the 5 molecules in Fig. 6. In addition, the individual time-dependent states are reported in Table III and if they fell outside of the Franck-Condon spectrum, their energies are marked by the red stick absorption on Fig. 6.

State . | TDDFT energy . | f
. | Assign. . | Sym. . | Expt. (eV) . |
---|---|---|---|---|---|

ZnTBP | |||||

1A_{1g} | Ground state | D_{4h} | |||

1E _{u} | 2.06 eV (601 nm) | 0.2550 | Q (π-π∗) | D_{4h} | 2.06^{a} |

2.02^{b} | |||||

2E _{u} | 3.56 eV (348 nm) | 1.475 | B (π-π∗) | D_{4h} | 3.18^{a} |

3.06^{b} | |||||

ZnTBMAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1B_{2} | 2.05 eV (605 nm) | 0.3684 | Q (π-π∗) | C_{2v} | |

1A_{1} | 2.13 eV (582 nm) | 0.2507 | Q (π-π∗) | C_{2v} | |

1A | 3.59 eV (345 nm) | 1.294 | B (π-π∗) | C_{2} | |

1B | 3.69 eV (336 nm) | 1.017 | B (π-π∗) | C_{2} | |

ZnTBcisDAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1A_{1} | 2.03 eV (612 nm) | 0.4184 | Q (π-π∗) | C_{2v} | |

1B_{2} | 2.17 eV (572 nm) | 0.2919 | Q (π-π∗) | C_{2v} | |

1B | 3.52 eV (352 nm) | 0.6420 | B (π-π∗) | C_{1} | |

1A | 3.84 eV (323 nm) | 0.7056 | B (π-π∗) | C_{1} | |

ZnTBtransDAP | |||||

1A _{g} | Ground state | D_{2h} | |||

1B_{3u} | 1.94 eV (640 nm) | 0.4742 | Q (π-π∗) | D_{2h} | |

1B_{2u} | 2.12 eV (584 nm) | 0.2722 | Q (π-π∗) | D_{2h} | |

1B | 3.49 eV (355 nm) | 0.8643 | B (π-π∗) | C_{2} | |

1A | 3.73 eV (333 nm) | 0.6223 | B (π-π∗) | C_{1} | |

ZnDPTBtransDAP | |||||

1A _{g} | Ground state | D_{2h} | |||

1A | 1.90 eV (650 nm) | 0.5196 | Q (π-π∗) | C_{1} | 1.85^{c} |

2A | 2.12 eV (586 nm) | 0.2408 | Q (π-π∗) | C_{1} | 1.96^{c} |

Q | 2.06^{c} | ||||

1B | 3.35 eV (370 nm) | 0.9145 | B (π-π∗) | C_{1} | 2.89^{c} |

2B | 3.56 eV (348 nm) | 1.1359 | B (π-π∗) | C_{1} | 3.82^{c} |

ZnTBTrAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1A | 1.94 eV (640 nm) | 0.4944 | Q (π-π∗) | C_{1} | |

1B_{2} | 2.07 eV (598 nm) | 0.4663 | Q (π-π∗) | C_{2v} | |

2A | 3.43 eV (361 nm) | 0.3024 | B (π-π∗) | C_{1} | |

3A | 3.86 eV (321 nm) | 0.2939 | B (π-π∗) | C_{1} | |

ZnMPTBTrAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1A | 1.93 eV (642 nm) | 0.5115 | Q (π-π∗) | C_{1} | 1.85^{c} |

2A | 2.12 eV (586 nm) | 0.3463 | Q (π-π∗) | C_{1} | 1.92^{c} |

Q | 2.04^{c} | ||||

3A | 3.40 eV (364 nm) | 0.4617 | B (π-π∗) | C_{1} | 2.81^{c} |

4A | 3.55 eV (349 nm) | 0.6613 | B (π-π∗) | C_{1} | 3.25^{c} |

State . | TDDFT energy . | f
. | Assign. . | Sym. . | Expt. (eV) . |
---|---|---|---|---|---|

ZnTBP | |||||

1A_{1g} | Ground state | D_{4h} | |||

1E _{u} | 2.06 eV (601 nm) | 0.2550 | Q (π-π∗) | D_{4h} | 2.06^{a} |

2.02^{b} | |||||

2E _{u} | 3.56 eV (348 nm) | 1.475 | B (π-π∗) | D_{4h} | 3.18^{a} |

3.06^{b} | |||||

ZnTBMAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1B_{2} | 2.05 eV (605 nm) | 0.3684 | Q (π-π∗) | C_{2v} | |

1A_{1} | 2.13 eV (582 nm) | 0.2507 | Q (π-π∗) | C_{2v} | |

1A | 3.59 eV (345 nm) | 1.294 | B (π-π∗) | C_{2} | |

1B | 3.69 eV (336 nm) | 1.017 | B (π-π∗) | C_{2} | |

ZnTBcisDAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1A_{1} | 2.03 eV (612 nm) | 0.4184 | Q (π-π∗) | C_{2v} | |

1B_{2} | 2.17 eV (572 nm) | 0.2919 | Q (π-π∗) | C_{2v} | |

1B | 3.52 eV (352 nm) | 0.6420 | B (π-π∗) | C_{1} | |

1A | 3.84 eV (323 nm) | 0.7056 | B (π-π∗) | C_{1} | |

ZnTBtransDAP | |||||

1A _{g} | Ground state | D_{2h} | |||

1B_{3u} | 1.94 eV (640 nm) | 0.4742 | Q (π-π∗) | D_{2h} | |

1B_{2u} | 2.12 eV (584 nm) | 0.2722 | Q (π-π∗) | D_{2h} | |

1B | 3.49 eV (355 nm) | 0.8643 | B (π-π∗) | C_{2} | |

1A | 3.73 eV (333 nm) | 0.6223 | B (π-π∗) | C_{1} | |

ZnDPTBtransDAP | |||||

1A _{g} | Ground state | D_{2h} | |||

1A | 1.90 eV (650 nm) | 0.5196 | Q (π-π∗) | C_{1} | 1.85^{c} |

2A | 2.12 eV (586 nm) | 0.2408 | Q (π-π∗) | C_{1} | 1.96^{c} |

Q | 2.06^{c} | ||||

1B | 3.35 eV (370 nm) | 0.9145 | B (π-π∗) | C_{1} | 2.89^{c} |

2B | 3.56 eV (348 nm) | 1.1359 | B (π-π∗) | C_{1} | 3.82^{c} |

ZnTBTrAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1A | 1.94 eV (640 nm) | 0.4944 | Q (π-π∗) | C_{1} | |

1B_{2} | 2.07 eV (598 nm) | 0.4663 | Q (π-π∗) | C_{2v} | |

2A | 3.43 eV (361 nm) | 0.3024 | B (π-π∗) | C_{1} | |

3A | 3.86 eV (321 nm) | 0.2939 | B (π-π∗) | C_{1} | |

ZnMPTBTrAP | |||||

1A_{1} | Ground state | C_{2v} | |||

1A | 1.93 eV (642 nm) | 0.5115 | Q (π-π∗) | C_{1} | 1.85^{c} |

2A | 2.12 eV (586 nm) | 0.3463 | Q (π-π∗) | C_{1} | 1.92^{c} |

Q | 2.04^{c} | ||||

3A | 3.40 eV (364 nm) | 0.4617 | B (π-π∗) | C_{1} | 2.81^{c} |

4A | 3.55 eV (349 nm) | 0.6613 | B (π-π∗) | C_{1} | 3.25^{c} |

The only spectrum from Fig. 6 that had comparable experimental data was ZnTBP, and as can be seen in Table III, the Q-band at the 0-0 transition was calculated at 2.06 eV for the degenerate state and was exactly the same as the jet gas phase. The B-band transition was again a degenerate state, but the energy was calculated at nearly 0.50 eV higher than experiment. While the Q-band results are close enough to experiment for ZnPc and ZnTBP, the rest of the azaporphyrin theoretical results would be questionable in the B-band. Since these energies can depend heavily on the Hartree-Fock exchange energy used as noted by Nemykin *et al.*, one could conclude that while the separated range method M11 was the best functional for the Q band, it may not be the best functional to look at the lower excited-state energies.^{26}

There was only one azaporphyrin analog where the symmetry was broken in the Q band, this was ZnTBTrAP in the first excited state at 1.9375 eV. As discussed in Sec. V A, the energy difference between the minimum and the barrier has to be investigated for reliability of the FC spectrum. The energy difference between the minimum and the barrier was 0.0207 eV. If the energy is doubled for an energy buffer above the minimum, this energy comes to 0.0414 eV. While this energy is a bit larger than what was seen in ZnPc, the second excited state is calculated at 2.07 eV, so there is again no reason to consider the FC result to be unreliable.

### C. Excited state theory compared to experiment for phenyl additions

ZnMPTBTrAP and ZnDPTB*trans*DAP were clearly two molecules of interest that came out of the molecular orbital levels as seen in Fig. 4 where the HOMO LUMO gaps were 2.18 eV and 2.16 eV, respectively. In addition, the absorption predictions in Fig. 7 show strong peaks in the Q band region for both of these molecules. The excited state calculations, as described in Secs. V A and V B, were carried out using the same method of TDDFT with the M11 functional. The excited state results can be seen in Table III, along with experimental data. Again, the excited state calculations over estimated the energies in the B Band, but the Q-band results were only 0.05 eV-0.10 eV different from experiment. Symmetry is reported in Table III for the ground state, but no symmetry was used in the calculations of the ground or excited states due to the high probability of the phenyl rings becoming distorted in the minimum energy optimization. However, even though no symmetry was used, the phenyl rings stayed orthogonal and were only slightly askew from an orthogonal position.

The Franck-Condon absorption spectra were also produced for these two molecules and overlaid with their non-phenylated sister molecules in Fig. 7. These two molecules were synthesized as described in Sec. III. The results in Fig. 7 show a better agreement with the experimental spectrum of ZnDPTB*trans*DAP as compared to theory, as well as a lower HOMO LUMO gap as compared to ZnMPTBTrAP, seen in Fig. 4. The results from ZnMPTBTrAP in Fig. 7 are more shifted with less definition in the Q band where there are two close sharp peaks of absorption in the experimental spectrum. In both cases as seen in Fig. 7, the two molecules have extremely similar shapes as compared to their non-phenylated pairs in the Q band region. Since the phenyl additions to these porphyrin analogs were orthogonal to the inner macrocycle, the fact that there was little difference in the spectra can be attributed to minimal interference of *π*-*π*^{∗} conjugation in the orbitals. This is also evident in Fig. 4, where there is little to no difference between the LUMO, HOMO, and HOMO-1 orbitals as compared with the phenylated pairs. The most promising result came from ZnDPTB*trans*DAP, where there was a strong absorption in the Q-band, and the HOMO LUMO gap reduced to 2.16 eV, lower than ZnPc.

## VI. CONCLUSION

The method of M11 worked well for the lowest energy excited state in all eight of the macrocyclic porphyrin analogs investigated. This separated range method did not do as well with the B band of the spectrum as the Q band, which could be due to the high percentage of the exchange energy at the long range. There would be a big advantage in the ability to tune the separated range for the center and the macrocycle. The 0-0 Q-band transition was predicted in vacuum as the exact same energy as gas phase for ZnPc. The prediction of a *C*_{4v} symmetry in the excited state Q band is a conclusion that can be made and fits with prior observed data. While the TDDFT and the Franck-Condon analysis over estimated the energy and oscillator strength for a few of the B bands, it worked well for every Q band. In general, although the time-dependent results were at times blue shifted, the Franck-Condon spectra that were produced were remarkably similar in shape to the few experimental spectra available in this research. The similarity in the absorption spectra with and without the phenyl rings is a key finding in this research. Saving time on these calculations by making a C–H bond assumption is a big advantage. Future work would include using other range-separated methods to see how the theoretical spectra might change. Needless to say it would be worth investigating a functional that would have a balance of results for both the Q band and the B band excited states as future work. Investigating absorption spectra that includes Herzberg-Teller terms would be included in future work. Making excited-state approximations for the phenyl additions that produced a buckled or saddle shape in the ground state would also be part of future work, in addition to looking at the excited states of these molecules in solvents. Overall, these molecules, especially ZnMPTBTrAP and ZnDPTB*trans*DAP, appear to be promising candidates for organic solar cells due to their similarity to ZnPc and a predicted reduction in recombination rates.

## Acknowledgments

Much thanks to Dr. Fernando Clemente for the tremendous help navigating the Franck-Condon calculations in Gaussian 09.

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant Nos. DGE 0802261 and DGE 1311230.

We would like to acknowledge and thank XSEDE and SDSC for a computational grant to work on Trestles, a dedicated XSEDE cluster designed by Appro and SDSC.