A real shift of 0.2 a.u. was used for the unsubstituted copper corrole and 0.4 a.u. was used for the triphenyl copper corrole. All CASSCF and CASPT2 energies were converged to 1.0 × 10 -8 a.u. The geometry optimization was considered converged when the maximum component of the gradient was 3.0 × 10 -4 a.u. per Bohr, the maximum displacement was 1.2 × 10 -3 Bohr, and the maximum change in energy was 1.0 × 10 -6 a.u. The SS-SR contraction scheme was used. The zeroth order Hamiltonian does not employ the so-called IPEA shift, and all computations were performed without imposing symmetry. To confirm that XMS-CASPT2 optimized geometries are minima, full harmonic vibrational analyses (FHVA) was performed by means of gradient finite difference for the S 0 , S 1 , and T 0 states of the bare copper corrole. Due to the large size of the triphenyl copper corrole, partial Hessian vibrational analysis (PHVA) was performed for the S 0 and T 0 states. The three phenyl groups were not included in the vibrational analysis. This allows for confirmation that the structure of the corrole itself is at a minimum. An in-house python code was then used to compute the molecular partition functions for the computation of standard state thermal contributions to Gibbs free energy at 298.15 K, employing the usual ideal gas, rigid-rotor, and quasi-harmonic oscillator approximations. The entropy and free energies are computed using the quasi-harmonic oscillator approximation proposed by Truhlar and Cramer in which all frequencies below 100 cm -1 are taken to be 100 cm -1 . For comparison, DFT gas phase optimizations are carried out employing the B3LYP, 81-84 PBE0, TPSSh, M06, BLYP, 83 PBE, TPSS, and M06-L 90 functionals using the def2-TZVP basis set on all atoms. The resolution of identity (RI) approximation was used for integral evaluation. The dispersion correction (D3) with Becke-Johnson damping was employed in all cases except M06 and M06-L. DFT geometries were confirmed as minima by harmonic vibrational analysis. Gibb's energy and enthalpy were calculated by replacing all normal mode below 100 cm -1 by 100 cm -1 . Time-dependent density functional theory (TDDFT) single-point calculations were performed on the TPSSh-D3 geometries using same method.