A Theoretical Study of the Occupied and Unoccupied Electronic Structure of High- and Intermediate-Spin Transition Metal Phthalocyaninato (Pc) Complexes: VPc, CrPc, MnPc, and FePc

The structural, electronic, and spectroscopic properties of high- and intermediate-spin transition metal phthalocyaninato complexes (MPc; M = V, Cr, Mn and Fe) have been theoretically investigated to look into the origin, symmetry and strength of the M–Pc bonding. DFT calculations coupled to the Ziegler’s extended transition state method and to an advanced charge density and bond order analysis allowed us to assess that the M–Pc bonding is dominated by σ interactions, with FePc having the strongest and most covalent M–Pc bond. According to experimental evidence, the lightest MPcs (VPc and CrPc) have a high-spin ground state (GS), while the MnPc and FePc GS spin is intermediate. Insights into the MPc unoccupied electronic structure have been gained by modelling M L2,3-edges X-ray absorption spectroscopy data from the literature through the exploitation of the current Density Functional Theory variant of the Restricted Open-Shell Configuration Interaction Singles (DFT/ROCIS) method. Besides the overall agreement between theory and experiment, the DFT/ROCIS results indicate that spectral features lying at the lowest excitation energies (EEs) are systematically generated by electronic states having the same GS spin multiplicity and involving M-based single electronic excitations; just as systematically, the L3-edge higher EE region of all the MPcs herein considered includes electronic states generated by metal-to-ligand-charge-transfer transitions involving the lowest-lying π* orbital (7eg) of the phthalocyaninato ligand.


Introduction
Phthalocyanines (H 2 Pc) share with porphyrins (H 2 P), everywhere present "as far as the living world is concerned" [1], the same four nitrogen-based coordinative pockets. Even though H 2 Pc and its metal complexes (MPc) are not present in Nature, they have been attracting great interdisciplinary interest because their technological potential spans over a wide range of applications [2,3]. Besides traditional appliances, such as dyestuffs for textiles and inks [3], MPcs are currently used as intrinsic semiconductors, chemical sensors, organic light-emitting diodes, organic photovoltaic cells, thin-film transistors, materials for nonlinear optics, spintronics and laser recording [4][5][6][7][8]. Moreover, bio-inspired oxygenbinding MPcs have been shown as viable substitutes for precious metals in catalysts for the oxygen reduction reaction (ORR) in low-temperature fuel cells [9][10][11][12]. As intimate an understanding as possible of the origin, symmetry, and strength of the M-Pc interaction is then mandatory to enhance the efficiency of new MPc-based devices. As such, X-ray absorption spectroscopy (XAS) is unanimously recognized as a valuable tool to probe, element-selectively, the empty electronic structure of M complexes, the M coordinative environment, as well as the nature and the strength of the M-ligand interaction [13][14][15][16]. recently synthesized in extreme conditions [35], the remaining MPcs herein considered all have relevant catalytic applications, including the: (i) CO and NO oxidation as well as ORR (CrPc) [49][50][51][52]; (ii) oxidation reactions in homogeneous and heterogeneous phase (MnPc) [9]; (iii) N-alkylation [53], C-H amination [54], C-C bond formation [55], synthesis of esters [56] and oximes [57] as well as reduction [58], oxidation [59,60], and radical reactions (FePc) [61]. Moreover, bio-inspired oxygen-binding Fe-macrocycles are particularly appealing as a consequence of their ability to reversibly bind O 2 , a crucial step in processes such as respiration, photosynthesis or the ORR catalysis [62]. Thus, there is no doubt that a comprehension as intimate as possible of the M-Pc interaction is of paramount importance to design new MPc-like species devoted to specific purposes.

Computational Details
The ground state (GS) electronic and structural properties of title molecules have been herein investigated by exploiting the Amsterdam Density Functional (ADF) package [63] within the assumption of an idealized D 4h symmetry [17] (see Figure 1), by running spin-unrestricted, nonrelativistic DFT calculations, with generalized gradient corrections self-consistently included through the Becke−Perdew formula [64,65], by adopting a tripleζ with a polarization function Slater-type basis set for all the atoms, and by freezing the M 1s-2p AOs and the 1s AO of N and C atoms throughout the calculations. MPc optimized Cartesian coordinates are reported in Tables S1-S4 of the Supplementary Materials. would then be utterly inadequate, while they can be properly handled by exploiting the module ROCIS of the ORCA program package [42].
In this contribution, a homogeneous modelling of the VPc, CrPc, MnPc, and FePc L2,3edges' features is presented and discussed with the ultimate goal of providing a first-principle rationale of differences characterizing the M-Pc interaction along the investigated series. As such, it can be useful to mention that, with the exception of VPc, recently synthesized in extreme conditions [35], the remaining MPcs herein considered all have relevant catalytic applications, including the: (i) CO and NO oxidation as well as ORR (CrPc) [49][50][51][52]; (ii) oxidation reactions in homogeneous and heterogeneous phase (MnPc) [9]; (iii) N-alkylation [53], C-H amination [54], C-C bond formation [55], synthesis of esters [56] and oximes [57] as well as reduction [58], oxidation [59,60], and radical reactions (FePc) [61]. Moreover, bio-inspired oxygen-binding Fe-macrocycles are particularly appealing as a consequence of their ability to reversibly bind O2, a crucial step in processes such as respiration, photosynthesis or the ORR catalysis [62]. Thus, there is no doubt that a comprehension as intimate as possible of the M-Pc interaction is of paramount importance to design new MPc-like species devoted to specific purposes.

Computational Details
The ground state (GS) electronic and structural properties of title molecules have been herein investigated by exploiting the Amsterdam Density Functional (ADF) package [63] within the assumption of an idealized D4h symmetry [17] (see Figure 1), by running spin-unrestricted, nonrelativistic DFT calculations, with generalized gradient corrections self-consistently included through the Becke−Perdew formula [64,65], by adopting a triple-ζ with a polarization function Slater-type basis set for all the atoms, and by freezing the M 1s-2p AOs and the 1s AO of N and C atoms throughout the calculations. MPc optimized Cartesian coordinates are reported in Tables S1-S4 of the Supplementary Materials. Insights into the origin, symmetry, and strength of the M-Pc interaction have been gained by combining the Nalewajski−Mrozek approach [66][67][68][69][70][71], well suited to estimating Insights into the origin, symmetry, and strength of the M-Pc interaction have been gained by combining the Nalewajski−Mrozek approach [66][67][68][69][70][71], well suited to estimating bond multiplicity indices ( NM I) in M complexes [66][67][68][69][70][71][72], with the Ziegler's extended transition state (ETS) method [73]. According to the ETS scheme, the MPc bonding energy (BE) may be written as: where ∆E es accounts for the pure electrostatic interaction, ∆E Pauli represents the destabilizing two-orbital−four-electrons interaction between the occupied orbitals of the interacting fragments (only atomic fragments have been herein considered), and ∆E orb corresponds to the stabilizing interaction between the occupied and empty orbitals of the atomic fragments. In passing, ∆E orb may be further decomposed into contributions due to the different irreducible representations (IRs) of the D 4h point group, according to The MPc L 2,3 -edges' XA spectra [33][34][35] have been modelled by evaluating EEs and corresponding oscillator strength distributions (f (EE)) for transitions with the M 2p-based MOs as initial spin orbitals (isos), by means of the DFT/ROCIS method [30], which includes SOC in a molecular Russell−Saunders fashion [25], by adopting the B3LYP exchangecorrelation (XC) functional [74] for VPc, CrPc and FePc and the M06 meta-GGA XC [75] for MnPc (vide infra), and by using the def2-TZVP(-f) basis set [76,77]. The combined use of DFT and configuration interaction requires a set of three semi-empirical parameters (c 1 = 0.18, c 2 = 0.20, and c 3 = 0.40), which have been calibrated by Roemelt and Neese [25] for a test set of M L 2,3 -edges. Throughout the M L 2,3 -edges modelling, the resolution of identity approximation has been used with the def-TZVP/J basis set [76,77]. Moreover, the zero-th order regular approximation has been adopted to treat the scalar relativistic effects [78]. Numerical integrations for DFT/ROCIS calculations have been carried out on a dense Lebedev grid (302 points) [79]. In addition, MPc-modelled spectra have been shifted by 10.9 (VPc), 9.8 (CrPc), 8.3 (MnPc) and 13.1 (FePc) eV to superimpose the highest intensity features of the simulated and experimental L 3 -edge, which does not suffer from the extra broadening and the distortion due to the Coster-Kronig Auger decay process [32,80]. This was needed because absolute theoretical EEs carry errors arising from DF deficiencies in the core region, one-particle-basis set restrictions and inadequacies in the modelling of spin-free relativistic effects [24].

MPc Occupied Electronic Structure
The aim of obtaining an understanding as intimate as possible of the origin, symmetry and strength of the M-Pc interaction may benefit from a preliminary, qualitative description of the MPc frontier orbitals simply based on symmetry arguments and overlap considerations. MPcs are united by the presence of the Pc 2− ligand whose electronic properties have been thoroughly described elsewhere [81]. Pc 2− frontier MOs may be split into σ and π sets. MPc symmetry adapted linear combinations (SALC) of C and N 2p σ (C and N 2p π ) are bases for the following D 4h IRs: a 1g , a 2g , b 1g , b 2g , e u (e g , a 1u , a 2u , b 1u , b 2u ); moreover, among π MOs, no a 1u SALC of N 2p π AOs, no b 1u SALC of N Py 2p π AOs and no b 2u SALC of N m 2p π AOs (see Figure 1) is present. In addition, the four N Py lone pairs pointing toward the centre of the coordinative pocket are bases for the a 1g , b 1g and e u IRs. In more detail, the two Pc 2− highest occupied MOs (HOMOs) correspond to the 15b 1g MO, σ in character and strongly localized on the N Py lone pairs pointing toward the centre of the coordinative pocket, and the 2a 1u π MO. In this regard, it is noteworthy that the D 4h a 1u IR is anti-symmetric with respect to the reflections through the σ h , σ v and σ d symmetry planes of the D 4h point group; the a 1u π MOs have then a node on symmetry planes passing through N Py and N m atoms (see Figure 1). As far as the D 4h Pc 2− lowest unoccupied MO (LUMO) is concerned, the 6e g VMO has a π character too.
The presence of the Pc 2− square planar ligand field lifts the five-fold degeneracy of the M 3d AOs, generating a 3d σ and a 3d π set. The former set includes the 21a 1g (z 2 -based) and the 16b 1g (x 2 -y 2 -based) MOs, while the latter takes in the 14b 2g (xy-based) and the 6e g (xz-and yz-based) ones [17]. Among them, the 21a 1g and the 14b 2g MOs are substantially M-N Py non-bonding, while the 6e g and 16b 1g MOs are M-N Py π and σ antibonding, respectively. Relative energy positions of the 3d σ /3d π spin up (↑)/spin down (↓) sets in VPc, CrPc, MnPc and FePc are displayed in Figure 2, together with those of selected Pc-based MOs.
The presence of the Pc 2− square planar ligand field lifts the five-fold degeneracy of the M 3d AOs, generating a 3d and a 3d set. The former set includes the 21a1g (z 2 -based) and the 16b1g (x 2 -y 2 -based) MOs, while the latter takes in the 14b2g (xy-based) and the 6eg (xz-and yz-based) ones [17]. Among them, the 21a1g and the 14b2g MOs are substantially M-N Py non-bonding, while the 6eg and 16b1g MOs are M-N Py  and  antibonding, respectively. Relative energy positions of the 3d/3d spin up ()/spin down () sets in VPc, CrPc, MnPc and FePc are displayed in Figure 2, together with those of selected Pc-based MOs. VPc ground state. Experimental [35] and theoretical [36,37,82] results disagree about the VPc ground state (GS) spin multiplicity. Carlotto et al. [36,37] recently proposed a 4 Eg high-spin (HS) GS ( 1 1 2 1 1 1 0 ; see Figure 2 and   VPc ground state. Experimental [35] and theoretical [36,37,82] results disagree about the VPc ground state (GS) spin multiplicity. Carlotto et al. [36,37] recently proposed a 4 E g high-spin (HS) GS (a 1 1g b 1 2g e 1 g b 0 1g ; see Figure 2 and Table 1; possible Jahn-Teller distortions [83] associated with orbitally degenerate GS or excited states have not been taken into account), while Eguchi et al. [35] presumed a 2 E g low-spin (LS) GS (a 0 1g b 2 2g e 1 g b 0 1g ). The 2 Eg state generated by the a 0 1g b 2 2g e 1 g b 0 1g configuration is 52.1 (35.5) kcal/mol less stable than the HS 4 E g (LS 2 B 1g ) one; moreover, the 4 B 1g / 4 A 2g states generated by the constrained a 0 1g b 1 2g e 2 g b 0 1g /a 1 1g b 0 2g e 2 g b 0 1g configurations are 1.8/1.7 kcal/mol less stable than the 4 Eg GS. Consistently with the VPc a 1 1g b 1 2g e 1 g b 0 1g GS configuration, the NM I V-N Py is quite large (0.64) [36]; as such, since NM I includes both covalent and ionic contributions, it is of some relevance to mention that the V Hirshfeld [84] charge (Q V ) amounts to 0.34. A closer look at the frontier VPc GS electronic structure indicates that all the V 3dbased singly occupied MOs (SOMOs) lie well above the ring-based, V-free, 2a 1u π doubly occupied MO (DOMO). The ionization energies (IEs) of VPc frontier MOs are not available in the literature; nevertheless, their values may be estimated by exploiting the Slater transition state (TS) method [85], which allows the evaluation of excitation/ionization energies " . . . by means of an artificial state that is halfway between the ground state of an atom or molecule and an excited state" [86]. Interestingly, the lowest VPc TSIE (6.17 eV) is associated with the ionization from the V 3d π -based 6eg ↑ SOMO rather than with the photoemission from the ring-based, V-free, 2a 1u π DOMO (6.59 eV). In this regard, it is of value to highlight that in his seminal paper devoted to the investigation of gas-phase photoelectron (PE) spectra of H 2 Pc and MPc (M = Mg, Fe, Co, Ni, Cu, and Zn), J. Berkowitz pointed out that " . . . the first ionization potential occurs at~6.4 eV, and it varies almost imperceptibly from sample to sample, including metal free and MgPc. Therefore, the conclusion seems inescapable that the first ionization potential corresponds to electron ejection from a ring orbital, and not a metal orbital" [87].
The comparison of the VPc bond lengths and bond angles with those optimized for the other MPc (see Table 2) [88] indicates that the structural perturbations induced by the presence of different M II ions in the Pc 2coordinative pocket are rather minute.
CrPc ground state. As already mentioned, CrPc has been attracting great interest as a catalyst for the CO and NO oxidation, as well as for the ORR [49][50][51][52], thus making particularly interesting the study of its electronic structure. Cr II has a 3d 4 configuration, which may generate three spin states with S = 0 (LS), S = 1 (intermediate spin; IS) and S = 2 (HS). Any attempt to optimize the LS state failed (NC, non-converged in Table 1), while the 3 E u IS state, associated with the a 1 1g b 1 2g e 2 g b 0 ↓ configuration, has been found less stable than the 5 B 1g HS one, generated by the a 1 1g b 1 2g e 2 g b 0 1g configuration, by 25.6 kcal/mol (IS and HS CrPc optimized structures are perfectly superimposable). Incidentally, the 3 E u IS state implied a non-Aufbau occupation accompanied by a pseudo reduction (oxidation) of the Cr II ion (macrocycle). As such, even though the 5 B 1g HS GS has been experimentally revealed [89,90] and theoretically predicted [36,82], the localization of VMOs is still controversial. Indeed, SIESTA [91] numerical experiments carried out by Arillo-Flores et al. [82] are consistent with the absence of " . . . metal contributions to HOMO and LUMO, they principally localize upon the inner ring.", which is certainly correct for the M-free 2a 1u HOMO, but wrong for the 6e g ↓ LUMO. In addition, differently from VPc, the ring-based, Cr-free, the 2a 1u π DOMO corresponds to the CrPc HOMO (see Figure 2). Analogously to VPc, the IEs of the CrPc frontier MOs are not available in the literature, but, differently from VPc, the lowest TSIE value (6.60 eV) is estimated for the ionization from the ring-based, Cr-free, 2a 1u π HOMO. Upon moving from VPc to CrPc, the GS frontier electronic configuration evolves from . The addition of an electron to the 3d π -based 6e g ↑ MO, M-N Py anti-bonding could then be invoked to rationalize the NM I M-N Py reduction from 0.64 (see above) to 0.43. Nevertheless, three things need to be kept in mind before drawing conclusions: (i) as already mentioned, NM I includes both covalent and ionic contributions; (ii) Q Cr (0.49) is larger than Q V (0.34); (iii) the 3d π -based 6e g ↑ MO is more anti-bonding in CrPc than in VPc (the localization % of the VPc 6e g ↑ MO on N Py is negligible; see Figure 3).  Upon moving from VPc to CrPc, the GS frontier electronic configuration evolves from  Analogously to VPc, no crystallographic data are available in the literature for CrPc; nevertheless, the tiny differences between the CrPc-and VPc-optimized structural parameters seem to indicate that a subtle balance between ionic and anti-bonding covalent contributions to the M-Pc interaction, both of them larger in CrPc than in VPc, takes place.
MnPc ground state. Likewise CrPc, LS, IS and HS states are possible; moreover, even though the optimized structural parameters corresponding to different spin states are very similar, the 4 Eg IS state associated to the 1 1 2 1 3 1 0 configuration (see Figure 2) has been found more stable than the 6 Eg HS and the 2 B2u LS ones by 16.5 and 14.9 kcal/mol, respectively. As such, it has to be noted that: i) the MnPc IS GS has been experimentally [94][95][96][97][98][99] and theoretically [36,[100][101][102][103][104][105] assessed; ii) the 6 Eg HS state is generated by the 1 1 2 1 2 1 0 7 1 configuration; iii) the 2 B2u LS state has the following occupation numbers:  Analogously to VPc, no crystallographic data are available in the literature for CrPc; nevertheless, the tiny differences between the CrPc-and VPc-optimized structural parameters seem to indicate that a subtle balance between ionic and anti-bonding covalent contributions to the M-Pc interaction, both of them larger in CrPc than in VPc, takes place.
MnPc ground state. Likewise CrPc, LS, IS and HS states are possible; moreover, even though the optimized structural parameters corresponding to different spin states are very similar, the 4 E g IS state associated to the a 1 1g b 1 2g e 3 g b 0 1g configuration (see Figure 2) has been found more stable than the 6 E g HS and the 2 B 2u LS ones by 16.5 and 14.9 kcal/mol, respectively. As such, it has to be noted that: (i) the MnPc IS GS has been experimentally [94][95][96][97][98][99] and theoretically [36,[100][101][102][103][104][105] assessed; (ii) the 6 E g HS state is generated by the a 1 1g b 1 2g e 2 g b 0 1g 7e 1 g configuration; (iii) the 2 B 2u LS state has the following occupation numbers: numbers of the Mn 3d-based 21a 1g , 14b 2g and 6e g MOs. Three different configurations may be considered: a 1 1g b 1 2g e 3 g b 0 1g ( 4 E g ), a 1 1g b 2 2g e 2 g b 0 1g ( 4 A 2g ) and a 2 1g b 1 2g e 2 g b 0 1g ( 4 B 1g ). In agreement with the theoretical results of Brumboiu et al. [105], the ADF outcomes rule out the 4 B 1g state because of its high energy; moreover, it is noteworthy that different experimental studies support either a 4 E g or a 4 A 2g GS. Specifically, XAS evidence [98], magnetic circular dichroism (MCD)/UV-Vis results [95], and XAS/MCD outcomes [99] favour a 4 E g GS [107], while magnetic susceptibility measurements have been rationalized within the assumption of a 4 A 2g GS [95,97], which has been attributed to intermolecular interactions in the MnPc crystal. In agreement with the literature [101], the ADF results herein reported estimate the 4 E g state to be more stable than the 4 A 2g by 7.3 kcal/mol. Before going on, it is of some relevance to point out that the MnPc a 1 1g b 1 2g e 3 g b 0 1g GS configuration implies the presence of a high-lying Mn 3d π -based SOMO ↓ well above the ring-based, Mn-free, 2a 1u π DOMO (see Figure 2). As such, no gas-phase photoemission results are available in the literature for MnPc; however, Grobosch et al. [108] were able to record the He(I) photoemission spectrum of an MnPc thin film deposited on polycrystalline Au. Interestingly, they assigned the lowest lying peak of the MnPc valence band photoemission spectrum to the ionization from the Mn 3d π -based 6e g ↓ MO, observing also that the ∆IE between the ring-based, Mn-free, π 2a 1u DOMO and the 3d π -based 6e g ↓ orbital is~0.5 eV [108,109]. In perfect agreement with these findings [108], the TSIEs [85] of the highest-lying 6e g ↓ and 2a 1u ↓ orbitals are 5.99 and 6.51 eV, respectively. Among the investigated molecules, MnPc is the lightest one for which structural parameters are available [88,92]. The data reported in Table 2 reveal that optimized bond lengths and bond angles fairly reproduce experimental evidence. In this context, the NM I Mn-N Py value (0.52), just in between the NM I Cr-N Py and the NM I V-N Py ones (0.43 and 0.64, respectively), seems to indicate that the M-N Py bond-weakening associated with the addition of a further electron to the 3d π -based 6e g MO is negligible (see Figure S1 of the Supplementary Materials). In addition, it has to be underlined that Q Mn (0.33) and Q V (0.34) are almost identical.
FePc ground state. Similarly to MnPc, experimental [34,[110][111][112][113][114] and theoretical [38,104,115] evidence indicates a FePc IS GS whose symmetry is, however, still debated. On the computational side, the IS GS symmetry, inextricably linked to the occupation numbers of frontier MOs, has been found to be extremely sensitive to the adopted XC functional and basis set. Carlotto et al. [38] have proposed a 3 E g IS GS, generated by a a 1 1g b 2 1g e 3 g b 0 2g configuration, on the basis of numerical experiments carried out by employing the ORCA program package [42], by using the hybrid B3LYP XC-functional [74], by adopting the def2-TZVP(-f) basis set [76,77] and the c 1 , c 2 , and c 3 semi-empirical parameters 0.21, 0.49, and 0.29 (hereafter, old set), respectively. In passing, the a 1 1g b 2 1g e 3 g b 0 2g GS configuration is a consequence of the orientation of the molecule in the xy plane. The x and y axes of the framework they adopted point toward the N m atoms rather than toward the N Py ones (see Figure 1). The four N Py lone pairs were then bases for the a 1g , b 2g and e u IRs rather than for the a 1g , b 1g and e u ones, and the Fe 3d-based VMO accounting for the Fe-N Py σ anti-bonding interaction corresponded to the 14b 2g level rather than to the 16b 1g one. All the possible configurations compatible with a FePc IS GS have been herein considered. As such, even though the ADF 3 B 2g (a 1 1g b 1 2g e 4 g b 0 1g ) and the 3 E g (a 1 1g b 2 2g e 3 g b 0 1g ) IS states are less stable than the 3 A 2g one (a 2 1g b 2 2g e 2 g b 0 1g ) by minute amounts (1.4 and 1.1 kcal/mol, respectively), it is of some relevance to underline that the ∆E χ orb (χ = a 1g , b 2g and e g ) of the IS states differ by up to~230 kcal/mol (∆E configuration. Besides the non-Aufbau filling of the corresponding electronic levels, the latter 3 E g state is significantly less stable than the 3 A 2g GS (12.3 kcal/mol). Additional numerical experiments have been carried out to estimate the FePc BE of the LS and HS states. As far as the former is concerned, it may imply either the a 0 1g b 2 2g e 4 g b 0 1g configuration or the a 2 1g b 0 2g e 4 g b 0 1g one. The LS a 2 1g b 0 2g e 4 g b 0 1g frozen configuration generates a non-Aufbau energy level filling, and a BE lower (~32 kcal/mol) than that of the 3 A 2g GS; any attempt to get a converged BE value for the LS a 2 1g b 0 2g e 4 g b 0 1g frozen configuration failed. Analogous considerations hold for the HS state. FePc is the lightest MPc for which gas-phase photoemission spectroscopy data have been recorded [87]. According to experimental evidence, the TSIE of the ring-based, Fe-free, π 2a 1u DOMO is the lowest one. Incidentally, its value (6.51 eV) is very close to the one (6.49 eV) estimated by Liao and Scheiner [101], even though a non-relativistic approach has been herein adopted.
According to the experiment in Ref. [93], the shortening of the M-N Py distance upon moving from MnPc to FePc is correctly reproduced (see Table 2). In this context, it is of relevance to mention that Q Fe = 0.13; i.e., the smallest value along the investigated series. Such a result, coupled to the NM I Fe-N Py = 0.64, ultimately indicates that the Fe-N Py interaction is the most covalent among the molecules we took into account (look at the MPc ∆E orb values reported in Table 3). In addition, the presence of a strong Fe-N Py covalent σ interaction is consistent with the high-energy position of the 16b 1g ↑ VMO (see Figure 2), which accounts for the σ anti-bonding interaction between the Fe 3d x 2 -y 2 AO and the b 1g SALC of the N Py lone pairs.
Electric dipole-allowed transitions imply that [17] where Γ GS , Γ µ , Γ ES and Γ Sym correspond to the IRs of the MPc electronic GS, the dipole moment operator (a 2u + e u ) [17], the electronic excited state (Γ ES = Γ iso ⊗ Γ GS ⊗ Γ fso ; iso and fso stand for initial and final spin orbitals, respectively), and the totally symmetric representation of the D 4h point group (a 1g ), respectively. Equation (3) may then evolve to which implies that, within the adopted approximation, which reduces the complete oneelectron excited configuration space (1h-1p space) to the subspace where only the M 2p core electrons (transforming as a 2u + e u ) [17] are excited, the allowed electric dipole transitions are (a 2u → a 1g ) ⊥ (5) (a 2u → e g ) || (6) (e u → e g ) ⊥ (e u → a 1g /a 2g /b 1g /b 2g ) || (8) where the ||/⊥ symbols stand for parallel/perpendicular to the molecular σ h plane (see Figure 1). VPc L 3 -edge. To date, only Eguchi et al. have succeeded in synthesizing, even though in extreme conditions, surface-supported mono-and multi-layers of VPc, whose angledependent linearly polarized XA spectra at the V L 2,3 -edges are reported in Ref. [35]. As such, a thorough analysis of the || and ⊥ V L 2,3 -edges components of the VPc and OVPc XA spectra has been recently reported by Carlotto et al. [36,37]. With specific reference to the VPc complex, the authors were able to conclusively assess its spin state by modelling corresponding XAS features for both LS (S = 1/2) and HS (S = 3/2) states. A brief description of their HS results [36] is herein included to favour the comparison with the modelled spectra of diverse MPcs.
It has been already mentioned that the one-electron excitation pattern describing the V II final states in D 4h symmetry should be dominated by states which may have (see Table  S5 of the Supporting Materials) either a spin multiplicity equal to (∆S = 0) or lower/higher (∆S = ±1) than the GS one [122]. As such, it has to be underlined that the ORCA B3LYP HS GS corresponds to the 4 A 2g state generated by the a 1 1g b 0 2g e 2 g b 0 1g electronic configuration. B3LYP/ROCIS outcomes [36,37] indicate that both the || f (EE) and the ⊥ f (EE) distributions of the L 3 -edge mainly arise from states having ∆S = 0, while ∆S = ± 1 contributions are negligible. Moreover, the lowest-lying ||/⊥ L 3 1 features (see Figure 4a) are due to states generated by V 2p → SOMOs single electronic excitations. Incidentally, SOMOs correspond to the V 3d-based 6e g and 21a 1g orbitals, while states associated to coupled-single electronic excitations [25,32], mainly involving V 2p → SOMOs and SOMOs → VMOs excitations, contribute to ||/⊥ L 3 2 . SOMOs naturally correspond to the 21a 1g and 6e g MOs, while the VMOs are the V 3d-based 14b 2g (~90%) and the π * Pc-based 7e g (~10%) orbitals. In passing, the MPc 7e g VMO corresponds to the lowest-lying Pc-based π * orbital and metal-to-ligandcharge-transfer (MLCT) transitions in diverse MPc L 3 -edge spectra (vide infra) involve this VMO. No contribution from V 2p → 16b 1g excitations is provided to states associated with the VPc L 3 -edge. The agreement between the experimental evidence and B3LYP/ROCIS outcomes has been documented in detail elsewhere [36,37]; here, it is sufficient to underline that both the relative positions and relative intensities of spectral features are well reproduced by ||/⊥ f (EE). Major disagreements between theory and experiment usually affect the L 2 region [31]; nonetheless, it deserves mentioning that the B3LYP/ROCIS || L 2 1 -|| L 3 2 HS ∆EE (6.4 eV, see Figure 4a) fairly reproduces the experimental value (6.9 eV). Any further comment about the VPc L 2 -edge is herein avoided. CrPc L3-edge. To date, the only CrPc f(EE) distribution for the 2p excitations is the one recorded by Koshino et al. by exploiting ISEELS [33]. Their spectrum includes both the L3and the L2-edge, but the coarse EE scale they adopted prevents the possibility of revealing the presence of possible structures associated to them; moreover, no information is provided by the authors about the L2-L3 ΔEE. Similarly to VPc, the one-electron excitation pattern describing the Cr final states in the CrPc D4h symmetry is dominated by states which may have (see Table S6 of the Supporting Material) a spin multiplicity either equal to (S = 0) or lower/higher (S = ±1) than the GS one. The HS CrPc  f(EE) and  f(EE) B3LYP/ROCIS distributions are superimposed upon the CrPc ISEEL spectrum in Figure  4b. In the L3-edge region, / f(EE) consists of an intense peak (L3 1 ) with an evident shoulder on its higher EEs (L3 2 , EE = 1.40 eV). Moreover, the B3LYP/ROCIS outcomes also indicate CrPc L 3 -edge. To date, the only CrPc f (EE) distribution for the 2p excitations is the one recorded by Koshino et al. by exploiting ISEELS [33]. Their spectrum includes both the L 3and the L 2 -edge, but the coarse EE scale they adopted prevents the possibility of revealing the presence of possible structures associated to them; moreover, no information is provided by the authors about the L 2 -L 3 ∆EE. Similarly to VPc, the one-electron excitation pattern describing the Cr final states in the CrPc D 4h symmetry is dominated by states which may have (see Table S6 of the Supporting Material) a spin multiplicity either equal to (∆S = 0) or lower/higher (∆S = ±1) than the GS one. The HS CrPc || f (EE) and ⊥ f (EE) B3LYP/ROCIS distributions are superimposed upon the CrPc ISEEL spectrum in Figure 4b. In the L 3 -edge region, ||/⊥ f (EE) consists of an intense peak (L 3 1 ) with an evident shoulder on its higher EEs (L 3 2 , ∆EE = 1.40 eV). Moreover, the B3LYP/ROCIS outcomes also indicate that both ∆S = 0 and ∆S = −1 states, both of them associated with single electronic excitations, contribute to the ||/⊥ f (EE) distributions. In more detail, ||/⊥ L 3 1 features are caused by ∆S = 0 states generated by Cr 2p → SOMOs transitions, with the whole set of the half-occupied Cr 3d-based orbitals (6e g , 21a 1g and 14b 2g SOMOs) as fsos. At variance with that, the ||/⊥ L 3 2 shoulders include contributions from both ∆S = 0 (64%) and ∆S = −1 (30%) states. Quintet states have the same origin as those associated with ||/⊥ L 3 1 , while the triplet ones are related to Cr-based 2p → 16b 1g and MLCT 2p → π * Pc-based transitions. As such, it is noteworthy that states associated with the Cr-based 2p → 16b 1g transition only contribute to || L 3 2 . Analogously to VPc, any comment about the CrPc L 2 -edge is herein avoided as it cannot be unambiguously determined by experiment [32]; nevertheless, we underline that the B3LYP/ROCIS L 2 -L 3 1 HS ∆EE (8.6 eV, see Figure 4b) fairly reproduces the experimental value (7.9 eV).
MnPc L 3 -edge. The MnPc L 2,3 -edges' XA spectrum recorded by Koshino et al. [33] suffers from the same issues already mentioned for CrPc. Otherwise, the experimental evidence reported by Kroll et al. [34] for MPc (25 ≤ Z ≤ 30) is much more informative as a consequence of the EE range they showed for each single MPc, thus allowing the detection of the fine structure eventually contributing to spectral features; moreover, both || and ⊥ polarized XA spectra are reported for each MPc. In the forthcoming discussion, ISEELS outcomes of Koshino et al. will no longer be considered as a reference for the L 2,3 -edges XAS modelling herein presented. At least three well-evident and closely spaced peaks (herein labelled || L 3 1 , || L 3 2 and || L 3 3 , and lying at~640.0,~641.0 and~643.0 eV, respectively) contribute to the MnPc || L 3 -edge spectrum (see Figure 4c). In addition, the || L 3 3 's higher EE side is characterized by the presence of an evident shoulder at~644.0 eV. The || → ⊥ light polarization switching is accompanied by a significant relative intensity reduction in spectral features having EEs in between 642 and 644 eV, thus reducing the ⊥ L 3 -edge spectrum to the ⊥ L 3 1 and ⊥ L 3 2 peaks with comparable intensity and lying at 640.5 and~641.5 eV, respectively (see Figure 4c).
The MnPc B3LYP/ROCIS ||/⊥ f (EE) distributions, not herein included, poorly reproduces the L 3 -edge XA spectrum in terms of number of peaks, relative energy positions and relative intensities. A few years ago, Carlotto et al. [123] tested the efficiencies of several diverse XC functionals (non-hybrid, hybrid and hybrid meta-GGA) in reproducing the L 2,3 -edges' absorption spectra of Mn complexes. The use of the hybrid M06 meta-GGA XC functional [75] was found to be decisive for a detailed assignment of the Mn(acac) 3 (acac = acetylacetonato) L 2,3 -edges XAS features. Now, despite the inability even of the M06 XC functional to reproduce in detail the complex structure of the MnPc XA spectrum, in particular the presence of the closely spaced || L 3 1 and || L 3 2 peaks, the IS-corresponding modelling (see Figure 4c) provides a satisfactory agreement between experiment and theory. In more detail, theoretical results indicate that the lowest-lying feature of the || f (EE) distribution ( || L 3 1 , representative of the || L 3 1 and || L 3 2 experimental peaks; see Figure 4c) has to be associated to the quartet electronic states (∆S = 0) generated by single electronic excitations having the whole Mn 2p set as isos and the low-lying Mn 3d-based 21a 1g ↓ , 14b 2g ↓ and 6e g ↓ VMOs as fsos. At variance with that, electronic states with ∆S = 0 (58%) and ∆S = 1 (24%) contribute to the intense || L 3 3 feature, while electronic states with ∆S = -1 negligibly contribute to the ||/⊥ L 3 patterns. Former states (∆S = 0) imply Mn 2p → SOMOs (21a 1g , 14b 2g and 6e g ) and SOMOs → 16b 1g /7e g coupled-single excitations [25,32], while the latter (∆S = 1) are mainly generated by Mn 2p → π * Pc-based single electronic excitations. The ∆S = 1 Mn 2p → π * Pc-based 4b 2u single electronic excitation violates selection rules stated by Eqns 5-8. As such, it has to be mentioned that i) ORCA calculations have to be run in C 1 symmetry and ii) the f value of the ∆S = 1 Mn 2p → π * Pc-based 4b 2u single electronic excitation is very low.
The satisfactory agreement between ||/⊥ f (EE) distributions and experimental evidence [34] prompts us to further detail the proposed assignment. Despite the fact that electronic states generating the ||/⊥ L 3 1 features of Figure 4c are associated with Mn-based ∆S = 0 2p → 3d single electronic excitations involving the whole 2p set and the low-lying Mn 3d-based ↓ VMOs, it sounds reasonable that the (1a 2u → 21a 1g ) ⊥ and (1e u → 6e g ) ⊥ transitions contribute to the L 3 lower EE region more than the (1a 2u → 6e g ) || , (1e u → 21a 1g ) || and (1e u → 14b 2g ) || ones, while the opposite is true when L 3 3 is considered. In fact, once again in agreement with the experimental results of Kroll et al. [34], the electronic states determined by ∆S = 0, 1 coupled-single and single ⊥ polarized electronic transitions provide to the L 3 higher EE region a contribution significantly lower than those || polarized.
DFT/ROCIS calculations fairly reproduce both the L 2 -L 3 ∆EE (~10 eV) and the corresponding relative intensities; nevertheless, any detailed assignment of the L 2 feature is herein avoided as it is not unambiguously determined by experimentation [32,80].
FePc L 3 -edge. FePc has been the object of a huge number of L 2,3 -edges XAS studies [33,34,38,120,[124][125][126][127]. Among them, Bartolomé et al. [125] investigated the Fe magnetic moment switching in the catalytic ORR of FePc adsorbed on Ag(110) by combining the results of X-ray linear polarized absorption spectroscopy with those of X-ray magnetic circular dichroism at the Fe L 2,3 -edges. A detailed analysis of the || and ⊥ Fe L 2,3 -edges components of the FePc and FePc(η 2 -O 2 ) XA spectra has been recently reported by Carlotto et al. [38] by adopting a computational set-up slightly different from that herein employed. The adopted set of c 1 , c 2 , and c 3 semi-empirical parameters corresponded to that herein labelled old set; moreover, the saturation of the final-state manifold was obtained by considering forty nonrelativistic roots per multiplicity.
The comparison of the FePc ||/⊥ f (EE) distributions (see Figure S2 of the Supplementary Materials) with the literature's theoretical results [38] proves an evident disagreement. To disentangle the effects induced by the adoption of a particular set of semi-empirical parameters from those generated by the number of nonrelativistic roots per multiplicity, the ||/⊥ f (EE) distributions have been again evaluated by using the c 1 , c 2 , and c 3 old set (see Figure 4d). Besides minor differences, most likely due to the higher number of states and expansion vectors herein adopted in the iterative solution of the CI equations, the ||/⊥ f (EE) distributions obtained by adopting c 1 = 0.21, c 2 = 0.49, and c 3 = 0.29 substantially mirror the literature's ones [38].
According to Carlotto et al. [38], weighty contributions to the ||/⊥ f (EE) distributions arise from states having a spin multiplicity either equal to (∆S = 0) or higher/lower (∆S = ±1) than the GS one. More specifically, the ||/⊥ L 3 1 features lying at the lowest excitation energy (see Figure 4d) are both associated to triplet states (∆S = 0), and are generated by single electronic excitations with the whole Fe 2p set as isos and the low-lying Fe 3d-based 21a 1g ↓ VMO and 6e g ↓ SOMO as fsos. Perfectly in agreement with the well-evidenced experimental dichroism [34,38], ⊥ L 3 1 is much more intense than || L 3 1 , thus indicating that electronic states associated with (a 2u → a 1g ) ⊥ /(e u → e g ) ⊥ excitations contribute much more than the (a 2u → e g ) || /(e u → a 1g ) || ones to the lower EE region of the L 3 -edge. Excitations with ∆S = 0, −1 comparably participate in states determining ⊥ L 3 2 (the theoretical setup herein adopted makes negligible ∆S = 1 contributions). As such, ∆S = 0 electronic excitations involve the Fe-based 6e g ↓ and 21a 1g ↓ VMO, while the ∆S = −1 ones have an MLCT character and imply Pc-based π * VMOs. Moving to the analysis of the || L 3 2 and || L 3 3 features, the comparison with the XAS evidence indicates that their excitation energies are slightly overestimated with respect to the || L 3 1 one; moreover, electronic states with ∆S = 0 (71%), 1 (19%) and −1 (7%) contribute to them. In even more detail, both single (44%) and coupled-single (56%) excitations contribute to the prevailing ∆S = 0 states. The former involves Fe-based (2p → 21a 1g ↓ /6e g ↓ , 2p → 16b 1g ) and MLCT (2p → 7e g ) transitions, while the latter are Fe-based (2p → SOMOs and SOMOs → 16b 1g ) excitations, determining the states contributing to the || L 3 2 and || L 3 3 lower excitation energy sides. The DFT/ROCIS calculations [38] fairly reproduce both the L 2 -L 3 ∆EE (~13 eV) and the corresponding relative intensities; nevertheless, any detailed assignment of the L 2 feature is herein avoided as it is not unambiguously determined by experimentation [32,80].

Conclusions
The occupied and empty states of HS VPc, CrPc, IS MnPc and FePc have been thoroughly investigated by exploiting the original/homogeneous theoretical results and experimental evidence form the literature. The use of the Hirshfeld charges [84] coupled with the Nalewajski−Mrozek [66][67][68][69][70][71] approach ultimately indicates that, among the investigated molecules, FePc is characterized by the strongest and most covalent M-Pc σ interaction. Even though Slater's transition state calculations ultimately confirm the Berkowitz hypothesis that, for FePc, " . . . the first ionization potential corresponds to electron ejection from a ring orbital, and not a metal orbital", the extension of the method to lighter MPcs reveals significant differences, the rationale of which lies with the relative energy position of the Pc 2− -based and M II 3d-based occupied/half-occupied MOs. Insights into the MPcs' virtual electronic structure have been gained by revisiting XAS data from the literature in light of DFT/ROCIS calculations. The higher EE side of the MPc L 3 -edge XA spectra systematically includes states associated with MLCT transitions, most of them involving the metal 2p → Pc 2− -based 7e g π * VMO excitations; moreover, the same EE region of all but one of (VPc) the L 3 -edge XA spectra is characterized by the presence of electronic states associated with M-based 2p → 16b 1g excitations. The agreement between theory and experiment is satisfactory, but it required a "tuning" of the modelling set-up in terms of XC functionals and/or c 1 , c 2 , and c 3 semiempirical parameters. As a final consideration, we underline that the theoretical outcomes obtained for the HS MPc (VPc and CrPc) are a true challenge for the experimental community called upon to confirm or deny them.   Table S5: Symmetries of a 2p 5 3d 4/6 system; Table S6: Symmetries of a 2p 5 3d 5 system; Table S7: Symmetries of a 2p 5 3d 7 system.