The Role of Self-Interaction Corrections , Vibrations , and Spin-Orbit in Determining the Ground Spin State in a Simple Heme

Without self-interaction corrections or the use of hybrid functionals, approximations to the density-functional theory (DFT) often favor intermediate spin systems over high-spin systems. In this paper, we apply the recently proposed Fermi–Löwdin-orbital self-interaction corrected density functional formalism to a simple tetra-coordinated Fe(II)-porphyrin molecule and show that the energetic orderings of the S = 1 and S = 2 spin states are changed qualitatively relative to the results of Generalized Gradient Approximation (developed by Perdew, Burke, and Ernzerhof, PBE-GGA) and Local Density Approximation (developed by Perdew and Wang, PW92-LDA). Because the energetics, associated with changes in total spin, are small, we have also calculated the second-order spin–orbit energies and the zero-point vibrational energies to determine whether such corrections could be important in metal-substituted porphins. Our results find that the size of the spin–orbit and vibrational corrections to the energy orderings are small compared to the changes due to the self-interaction correction. Spin dependencies in the Infrared (IR)/Raman spectra and the zero-field splittings are provided as a possible means for identifying the spin in porphyrins containing Fe(II).


Introduction
Because transition-metal atoms can exist in multiple spin and charge states and because of their large naturally occurring terrestrial abundance, nature has chosen such systems to form the basis of many life-maintaining molecular systems.For example, Fe-and Cu-porphyrin-like macrocycles carry oxygen for mammals and crustaceans, respectively, and Photosystem II, the oxygen evolving complex, is composed of spin-carrying Mn centers.For natural molecular systems, the chemical functions are impacted by varying spin states and/or charge states of the transition-metal system, and this is particularly important for the mammalian oxygen-transporting mechanisms that are accomplished by Heme-like molecules [1][2][3][4][5][6][7][8].A seemingly important characteristic of these molecules is that the multiple spin-and charge-states have similar degrees of stability/instability at room temperature leading to flexible transitions between different states.While the solvent environment is especially important for transitions involving changes in metal-center charge states, the on-molecule attributes allowing for low-energy spin-dependent transitions spring from the near energetic degeneracies of the transition-metal ions as a function of local moment for a specific charge state.Some evidence of energetic ambivalence in such systems may be deduced by noting the similarities and minor deviations in structure associated with the Mn centers that are predominantly undercoordinated in photosystem II but predominantly six-fold coordinated in the Mn 12 -Acetate molecular magnet [9,10].Synthesis of molecular magnets with six-fold coordination, accomplished by Lis et al. [11] and Christou et al. [12], suggests an energetic preference associated with six-fold coordinated centers.More recent work reinforces the suggestion that the six-fold coordination of transition metal centers [13,14] is favored over the lower-coordination systems.
With a few notable exceptions, such as the biologically relevant Fe 8 (TACN) molecule [15], where TACN is the organic ligand 1,4,7-triazacyclononane, and some Ni-cubane structures [16,17], the generalized-gradient approximation has provided very accurate predictions on the spin moments, spin-ordering, electronic structures, and magnetic anisotropy energies for molecular magnets composed of 4-12 six-fold coordinated high-spin transition metal centers that are held tightly in place by the spectator ligands.A combination of the charge states and shapes of the ligands create environments around the transition metal ions that mandate very specific d-state fillings and spin-states for each of the metal ions.For the aforementioned exceptions, Park et al. and Chao et al. have shown that Hubbard-U approximations can correct the local electronic structure of the transition metal ions [16,17].For the case of the Fe 8 (TACN) molecule, the early evidence for the failure to achieve a qualitatively correct electronic structure has been attributed to accidental degeneracies at the Fermi-level that lead to qualitatively incorrect electronic structures [15].More recently, Kulik et al. [18,19] have studied a wide class of single-center Fe compounds within standard and hybrid density-functional theory (DFT) methods as a function of the percentage of exact exchange and have shown that there can be strong variations in the energetically preferred iron spin state as a function of functional.In addition, Johnson et al. [20] have demonstrated similar uncertainties in moment size for single-iron molecular systems.For recent work on the role of exact exchange in hybrid functionals on the same molecule as discussed here, please see [21].As discussed below, Radon [21] has shown that exact exchange enhances the stability of high-spin transition-metal molecules and that both direct and self-consistent effects are important, which is consistent with the results presented below.
It is clear that an efficient first-principles methodology, with no adjustable parameters, is desired to more fully understand the energetics as a function of spin of under-coordinated metal ions.It is generally accepted that moment misorderings in transition-metal systems, when they exist, arise due to the self-interaction error [22][23][24][25] that is present in most energy functionals.The purpose of this paper is to unambiguously demonstrate that the self-interaction error impacts computational abilities to predict local spin moments and to demonstrate that a recently proposed formulation of the self-interaction-correction offers a path toward constraining future functionals to be self-interaction free and unitarily invariant.
The relative flexibility of the spin-state in the Fe(II) substituted porphyrin, hereafter referred to as Fe(II)-porphyrin, is evident from a lot of literature that shows uncertainties in the spin-ground state in these systems [1][2][3][4][5][6][7][8]26].Heme, an Fe(II) centered porphyrin group, is responsible for the oxygen transport.The Fe(II) ion, without O 2 binding, is about 0.6 Å above the center of the porphyrin ring.O 2 binding changes the heme's electronic structure, which shortens the Fe-N p bonds by ∼0.1 Å and Fe(II) moves into the heme plane.In vivo, the Fe(II) center in the porphyrin molecule is five-fold coordinated with the nitrogen on the apex of an imidazole acting as the fifth coordinating ligand.In the de-oxygenated state, the ground state of the complex is S = 2.The weak bonding between the imidazole and the porphyrin underscores the need to be able to computationally predict the spin states of the de-imidazolated porphyrin molecule also.If, as predicted by many density functionals, the de-imidazolated porphyrin has an S = 1 ground-state, the effective bonding between the porphyrin and the imidazole can be viewed as stronger since removal of the imidazole would then be spin forbidden.We concentrate on the spin states of the Fe(II)-porphyrin structure in this work.This is the simplest possible heme system to investigate within the Fermi-Löwdin-orbital self-interaction-corrected (FLOSIC) method.
In Section 2, the computational method for implementing and applying the Fermi-Löwdin self-interaction-correction to density functional theory, within the Naval Research Laboratory Molecular Orbital Library (NRLMOL), is reviewed.In Section 3, we present results on several different geometries and spin configurations of a simple mimic of Fe(II)-porphyrin molecule (FeN 4 C 20 H 12 ) using the PBE-GGA and PW92-LDA energy functionals.We show that inclusion of self-interaction corrections provides qualitative changes in the energetic ordering as a function of Fe(II) spin.In Section 4, we address additional interactions that could potentially impact the spin-dependent energetics of high-moment systems.We find that the overall weakening of bonds in the high-spin system decreases the zero-point energy of these systems relative to their lower-spin siblings and therefore slightly increases the stability of the high-spin Fe(II)-porphyrin system.We show further that the inclusion of spin-orbit energy ever-so-slightly increases the stability of the higher spin Fe(II)-porphyrin system relative to the larger-gap lower-spin systems.For the system studied here, the self-interaction correction is the major contributor to the qualitative change in stability as a function of the moment.Furthermore, since the SIC favors an Fe(II)-porphyrin with a permanent dipole, additional stabilization from solvents can be expected.

Theoretical and Computational Details
The calculations presented here have been performed with a modified version of the massively parallel NRLMOL density-functional suite of codes [27].This code uses very large basis sets [28] and a variational integration mesh [29] to perform numerically precise calculations on molecules composed of non-relativistic atoms with energies that are converged with respect to basis sets.Unless otherwise stated in the discussion, the molecular geometries are generally optimized until the Hellman-Feynman-Pulay forces [30] are smaller than 0.001 Hartree/Bohr, and Kohn-Sham orbitals used in the FLOSIC calculations have been self-consistently determined from the PBE-GGA functional.The vibrational results have followed the methods of Ref. [31].All geometries discussed in this work are found to have positive-definite hessian matrices.The corrections due to the spin-orbit interaction have been determined using the second-order post-SCF approaches that has been discussed in Refs.[9,10] to determine magnetic anisotropy energies.The one difference between the results reported there, which concentrate on the dependence of the spin-orbit energy as a function of spin orientation and the dependence of spin-orbit energy as a function of moment, is that one needs to include all occupied orbitals to determine the total spin-orbit energy that is needed here.
Generally speaking the best approach is to fully optimize a geometry but as discussed by Hahn et al., there may be cases where density-functional approximations yield a qualitatively incorrect geometry [32].However, it is generally posited that the S = 2 geometry for the Fe(II)-porphyrin is bent rather than planar.Collman et al. suggested that the intermediate-spin iron(II) atom should be precisely centered among the four complexed nitrogen atoms with an Fe-N bond distance smaller than 2.01 Å which is in accord with the geometries found here that have an Fe-N bond is 1.991 Å. Collman et al. state further that the high-spin iron(II) atom should lie substantially out of plane from the nitrogen atoms (probably by more than 0.45 Å) with an Fe-N bond distance larger than 2.07 Å.As such, we have created a high-spin iron geometry that is optimized to have the iron displaced by 0.38 Å above the plane and an Fe-N bond of 2.08 Å.Since most local and gradient-corrected functionals, as well as the ones used here, predict that there is no Jahn-Teller distortion or bending of the S = 2 geometry, we have also performed calculations on a bent S = 2 geometry to address the possibility that the self-interaction corrected functional will yield geometries that conform to experimental observations.

Fermi-Löwdin Orbitals for Unitarily Invariant SIC
The need for the use of self-interaction-corrections for systems containing transition-metal ions was first suggested by Lindgren [22] in applications to the Cu atoms in 1971.In 1981, Perdew and Zunger formalized the expression for the self-interaction correction, at the level of the total-energy expression, so that it could be applied to functionals designed to include correlation and exchange effects [23].Shortly thereafter, the Wisconsin SIC group [33][34][35] suggested the use of localization-equations as a variational means for addressing the unitary non-invariance that existed in the original formulation.This approach has been the primary means for variational self-interaction corrections as discussed in recent reviews [36][37][38][39][40][41][42][43].While the original version of SIC [33][34][35] only considered real orbitals, several researchers [40,44,45] have noted that there is no a priori reason to expect that complex orbitals could not lead to lower SIC energies and have derived additional formulae for their minimization.Pederson and Perdew noted in Ref. [46] that real localized orbitals that minimize the SIC energy provide stationary energies with respect to infinitesimal complex unitary transformations but again did not guarantee that the the real solutions could not be a saddle point or local maximum.It has been further pointed out that, if the SIC energy for localized orbitals is positive, delocalized orbitals, with vanishing SIC, would minimize the energy and it has been further noted that such behavior leads to a dilemma in the separated atom limit.
Recently, Pederson, Ruzsinszky and Perdew have reformulated the self-interaction-corrected density-functional expression so that it is explicitly unitarily invariant and size extensive [47].In this formulation, the SIC-energy depends on a set of Löwdin-orthogonalized Fermi-Orbitals (FLO) that are constructed from the spin-density matrix and Fermi-orbital descriptors (FODs), which often look like semiclassical electronic positions [32,[47][48][49].
In other words, the Fermi Orbital (FO) [50,51] depends parametrically on the variational position descriptor, a iσ , and defined according to: Because the FO accounts for all the spin density at a given point in space, it is expected to be a rather localized function.A picture of the locations of the Fermi-orbital descriptors may be found in Figure 1.The Löwdin symmetric reorthonormalization procedure guarantees that each localized orbital {φ iσ } is similar to its parent Fermi Orbital in a least-squares sense, but each localized orbital now depends parametrically on all the Fermi-orbital descriptors and on the spin-density matrix, or, equivalently, on the Kohn-Sham orbitals.Early indications [32,43,[47][48][49][52][53][54] are that the FLO-SIC improves orbital energies and atomization energies in a wide class of systems.
In general, the SIC formulation [23] proceeds by attaching an orbital-by-orbital self-interaction correction to whatever approximation to the universal functional is being employed.The SIC expression is given by: The terms U[ i,σ ] and E approx xc [ i,σ , 0] are the exact self-coulomb and approximate self exchange-correlation energies, respectively.In the above equation, the FLO {φ iσ } are used to define orbital densities according to: The derivatives of the energy with respect to the FOD can be derived analytically from the following expressions: with A complete derivation for the above terms and a discussion of the complexity can be found in Refs.[43,48].Given the derivatives, it is possible to optimize the SIC orbitals using conjugate-gradient or other gradient-based minimization procedures.

Method for Probing the SIC Energy Landscape
Obtaining starting guesses for the Fermi Orbital descriptors is initially difficult but transferable from one system to the next.Included in our optimization scheme here have been random moves of the FOD particularly for those associated with the macrocycle.To initially screen such moves, we use both total energy and the lowest eigenvalue of the Fermi Orbital overlap matrix as measures of probable success.
The self-interaction-correction used here as well as other formulations of self-interaction corrections lead to qualitative changes in the orbital energies.Because the self-interaction correction is generally negative, a self-interaction-corrected functional significantly discourages the possibility of finding electronic structures with fractionally occupied orbitals at the Fermi level.This in turn tends to pull all the fully occupied orbitals downward in energy.In both the original version of SIC and this version, this stabilization of the occupied orbitals improves the first and second derivatives of energy with respect to occupation numbers.The first derivatives are pushed downward leading to Janak's theorem, which is more similar to Koopman's theorem.The second derivatives of energy decrease significantly, which leads to a favoring of the much needed derivative discontinuity at integer occupations.Furthermore, by attenuating the possibility of fractionally occupied solutions, the inclusion of SIC strengthens the system's propensity toward Jahn-Teller distortions because it is no longer energetically possible to accept fractionally occupied solutions at the high-symmetry configurations.Jahn-Teller distortion, which explains the sacrifice of symmetric configuration of a molecule when orbital electronic degeneracy exists are known to be important in understanding the structure of porphyrins in the heme [55].

Spin-Orbit Coupling
In Refs.[9,10], equations are derived for calculating the magnetic anisotropy.The methods discussed in these papers have proven to provide very accurate means for determining the magnetic-anisotropy energies, also referred to as zero-field splittings, in molecular magnets.A primary reason for success is that the spin-carrying metal centers in such molecules are six-fold coordinated by ligands that strongly influence the d-shell filling in such systems.For specific applications to Fe(II) centers, the methods have been very successful in predicting the zero-field splitting of the Fe 4 molecular magnet [56][57][58] but less successful in obtaining quantitatively accurate zero-field splittings in the larger Fe 8 -TACN molecule [15] due to a qualitatively incorrect electronic structure.In this work, we have included the effects of spin-orbit interactions for several reasons.
First, the spin-orbit interaction contributes to the total energy of the system.While such effects are generally second order, some of the systems discussed here have vanishing highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) gaps or nearly vanishing HOMO/LUMO gaps.In such systems, a permanent orbital moment arises and this leads to a much stronger first-order contribution to the spin-orbit energy.Our results indicate that the permanent orbital moment does not lead to a strong enough spin-orbit energy to prevent a Jahn-Teller distortion.Second, the resulting magnetic anisotropy can be directly compared to experiment and used as an additional means for determining whether the spin state of the iron complex is S = 0, S = 1 or S = 2.These results are discussed in the individual sections and are tabulated in Table 2.
For this calculation of the spin-orbit contributions to the total energies and the zero-point splittings, the procedure developed in Ref. [9] is followed.For a collection of Kohn-Sham orbitals expressed as a product of a spatial orbital and spinor, ψ iσ ( r)χ σ , the second-order energy shift ∆ due to the spin-orbit interaction is determined from second-order perturbation theory according to: where σ and σ are spin indices, and iσ ( jσ ) are the eigenvalues of the occupied and unoccupied states, respectively.While the second-order expression derived in Ref. [9] is fine for the lowest-energy configurations here, for systems with permanent orbital moments, or partially occupied states at the Fermi level, the spin-orbit interaction needs to be determined from exact diagonalization methods.

Fe(II)-Porphyrin Structures vs. Spin, Orbital Moment and Spin Ordering
A characteristic of probable importance for the function of the Fe(II)-porphyrin is that there are several low lying spin states and in some cases there are multiple states of a given spin with different fillings of the d states.The possibility of multiple d-fillings was determined from the lower-symmetry calculations, discussed in the section on SIC.Based on that observation, we decided to perform calculations on the high-symmetry perfectly flat structures, which provide good reference calculations.Unless stated otherwise, all calculations discussed in this section are based upon the PBE-GGA energy functional.The full symmetry of the flat porphyrin can be generated by (x, y, z) → (y, x, z), (x, y, z) → (x, y, −z) and (x, y, z) → (x, −y, z).In Table 1, the equilibrium positions are provided for one S = 0 geometry, two S = 1 geometries and one S = 2 geometry.
Table 1.Spin-dependent structural and energetic data for flat D 4d Fe(II)-porphyrin structures are obtained with PBE-GGA.First row (left to right) shows the relative energies (E in eV), HOMO/LUMO gaps (∆λ in eV), the zero-field splitting (ZFS in Kelvin; EP/EA designates easy plane/axis), angular momentum (L) and spin-orbit energy (E LS in Kelvin).Jahn-Teller unstable cases are designated by JT n (n the number of degenerate states).Geometrical information follows (Å).Then, occupation numbers (superscript), symmetry label and orbital energies (eV) for frontier orbitals are included.All decimal places are included for data-management purposes.Fe atom is at the origin in all cases.The highest possible symmetry for the molecule is a perfect flat geometry with four-fold rotations about the z-axis.When the molecule is oriented with the nitrogen atoms at (A,0,0), (0,A,0), (−A,0,0) and (0,−A,0), (A = 1.98 Å), the Fe d-states are split into three one-fold representations (d xy , d z 2 , d x 2 −y 2 ) and one two-fold representation (d xz ,d yz ).The lobes of d x 2 −y 2 orbital are arranged to collide with the nitrogen atoms.As such, it is expected that both the kinetic energy of this orbital and the Coulomb repulsion between an electron in this orbital and electrons on the nitrogen atoms will make the d x 2 −y 2 orbital the least favorable orbital energetically, especially if the nitrogen has lone pairs pointing along the N-Fe axes.A picture of the d x 2 −y 2 in Figure 2 illustrates the penetration of the d-orbitals into the space that belongs to the nitrogen electrons.In contrast, the d xy and d z 2 orbitals arrange their lobes so that they avoid the nitrogens.From the point of view of coordination, the two-fold d xz /d yz states are less intrusive than the d x 2 −y 2 and might even form π bonds with nitrogen valence electrons.In fact, for this case, an unoccupied delocalized state, primarily associated with the pentagonal rings, is found to be lower than the localized d x 2 −y 2 state.For this configuration, we believe that the nearly degenerate d z 2 and d xz /d yz are higher in energy because the partial overlap of these densities increases the intra-shell Coulomb energy.Geometrical parameters and the values of orbital energy for this structure may be found in the fifth panel of Table 1.This is the highest energy geometry within both PBE-GGA and PW92-LDA.The exact degeneracy of the half-occupied d xz /d yz states, coupled with the near degeneracy with the fully occupied d 2 z states, allows the molecule to lower its energy through Jahn-Teller distortions and changes in total spin.A picture of the locations of the Fermi-orbital descriptors may be found in Figure 1.
Two S = 1 geometries (µ B = 2): From the spin unpolarized state, an electron configuration with four electrons occupying the majority spin d xy , d z 2 and d xz /d yz states may be created in two ways.In one case, found to be the lowest in energy, the electron can be removed from the minority spin d xz /d yz channels, which leads to a closed shell electronic configuration that is Jahn-Teller stable.In the second case, the electron can be removed from the d z 2 state, which leads to a higher energy albeit JT-unstable structure.The data may be found in the first and second panels of Table 1.
An Orbitally-Ordered Structure (µ B = 0): Starting from the spin unpolarized geometry, a lower energy state with no net moment may be found by transferring an electron from the d 2  z state and placing it in the hole of the d xz /d yz manifold of the same spin.In doing so, the Coulomb energy of the remaining d z 2 orbital decreases since there is no-longer the Coulomb energy between the now unpaired d 2 z state.This state may also be viewed as the partner of the open-shell S = 1 state, as it can also be generated by flipping the spin of one of the majority spin d xy /d yz states.This electronic configuration should be viewed as a mixture of open-shell singlet and triplet states.The data are listed in the third panel of Table 1.
The lowest S = 2 State: Analysis of this state suggests that it is a state that will be strongly affected by inclusion of the self-interaction correction.For the flat geometry, the lowest-energy S = 2 state has a closed majority spin d shell and has one electron in the d z 2 minority spin state.This state can be derived from the lowest energy S = 1 state by promoting an electron from the minority spin d xy state to the unoccupied majority spin d x 2 −y 2 state, which as mentioned above, is the LUMO+1 state of the lowest S = 1 state.When the iron atom relaxes out of plane, this state is able to lower its kinetic energy since its lobes are no longer directly penetrating the regions of space occupied by the nitrogen atoms.However, within PBE-GGA and PW92, we find that the flat pure S = 2 state is vibrationally stable and higher energy than the flat S = 1 states by 0.68 eV and 0.54 eV respectively.This reference configuration should not be as susceptible to multi-configurational corrections as the S = 1 states.It can be viewed as a pure S = 2 M s = ±2 state for either in-plane or out-of-plane configurations.This state is, however, expected to be strongly impacted by the self-interaction error.Its excess self-Coulomb repulsion tends to de-localize the d x 2 −y 2 state and the de-localization raises kinetic energy of the state and increases the mutual Coulomb interaction due to the fact that, spatially, the d x 2 −y 2 density resides in the regions close to the nitrogen electrons.

The S = 2 Porphyrin Structures: Flat vs. Bent
In addition to the flat geometry, discussed above, we have considered the configuration that the de-oxygenated geometry is found in vivo.We have optimized the geometry of the Fe(II)-porphyrin molecule, constrained to be in the high-spin (S = 2) state, using the PBE-GGA with the symmetry point group of order 4 that is generated with the symmetry operations (x → −x) and (y → −y).We have verified that PBE-GGA leads to a perfectly flat vibrationally stable structure and we have not found an out-of-plane S = 2 structure that has a lower energy within the PBE-GGA or PW92-LDA functionals.This finding is in accordance with the work of Kozlowski et al. [55].However, by analyzing the experimental geometries and the calculated vibrational modes, we have constructed an S = 2 out-of-plane porphyrin structure that is similar to experimental reports and is only 0.07 eV higher in energy than the perfectly flat D 4d structure.This out-of-plane geometry has an Fe ion that is 0.4136 Å above the planar average.There is a very small but non-zero gap (0.19 eV) between HOMO and LUMO indicating that this is an electronically stable configuration.The energy difference (0.025 eV) between the out-of-plane and perfectly-flat S = 2 geometries is smaller within the PW92-LDA functional.While the energy difference between the flat and bent structures are small, the orbital ordering in the minority spin channel changes.Rather than the occupation of the minority d z 2 orbital, found in the flat D 4d structure, the d yz orbital is pulled to lower energy through a series of mechanisms.First, simply promoting an electron from a fully occupied d z 2 state to a fully occupied d yz allows for a Jahn-Teller distortion.The fully occupied d xz orbital of opposite spin then hybridizes slightly with pentagonal rings on the porphyrins apparently to counteract increases in mutual Coulomb energy associated with the change in minority spin occupations.The three orbitals that are most important for the S = 2 structure are shown in Figures 2-4.One computational fingerprint associated with emptying of the localized minority spin d z 2 orbital and filling of the delocalized d yz orbital is that the total amount of electronic charge near the iron atom decreases.This slight change of charge state, coupled with clear rearrangement of orbital densities in the region that would be occupied by an oxygen molecule clearly argues for the need to account for these details.Here, we note that the energies of the two majority spin orbitals are directly affected by the inclusion of self-interaction corrections and further localization, due to full self-consistent treatments, is anticipated to be very important for the energetics of this system.With the PBE-GGA energy functional, we find that both S = 2 geometries are higher than two electronically distinct, but nearly degenerate, S = 1 structures (discussed below).Within PW92-LDA, the energy of the S = 2 structures are significantly higher than the the ground-state S = 1 structure.When the self-interaction correction is added, we find that the S = 2 geometry is lower, by 0.055 eV, in energy than either of the S = 1 geometries.When zero-point energies and spin-orbit energies are included, the stabilization of the S = 2 bent geometry is increased by an additional 0.05 eV.The 0.1 eV energy difference between the S = 2 and S = 1 geometry, for the self-interaction corrected PW92 result, may be susceptible to changes due to many small details of the calculation.However, the 1.3 eV stabilization of the S = 2 state relative to the S = 1 state suggests that the self-interaction correction is important in this molecule (See Table 2).
Table 2. Relative electronic energies as a function of functionals and spin for the lowest energy geometry of each spin.Also included are the zero-point vibrational energies and spin-orbit energies as calculated within PBE-GGA.All energies are in eV.Within FLOSIC-PW92, the S = 2 bent geometry is the lowest energy configuration found so far with an energy that is 0.055 eV lower in energy than the lowest-energy S = 1 geometry.Including zero-point effects further stabilizes the S = 2 geometry by approximately 0.050 eV.Inclusion of spin-orbit energy does not impact the energetic ordering.PW92-LDA and PBE-GGA find that the S = 1 state is higher in agreement with other works.In a recent paper [21], Radon has revisited the role of exact exchange in DFT spin-state energetics of transition-metal complexes including the FePy molecule discussed here.In this work, the author notes that "by increasing the amount of exact exchange, higher spin states are strongly stabilized relative to lower spin states".However, it is argued further that this effect is not primarily due to direct onsite lowering of the d-states but is substantially dependent on ligand-metal bonding, which suggests that there is also an indirect effect stated to be more important by Radon, which depends upon the relative ligand and metal frontier-orbital energies and the reduction of overlap which comes from direct accounting of orbital self-interaction on all sites.Such effects would results from self-consistent changes in the electronic wavefunctions.This observation is consistent with that of Chao et al. [17] who found that inclusion of Hubbard U corrections for both Ni 3d and O 2p state are important in similar systems.Radon finds that the S = 2 (quintet) state is higher in energy by 0.8 eV if full Slater exchange is used in the B3LYP functional but lower in energy by 0.8 eV when full exact exchange is used instead.For hybrid exchange calculations between these two endpoints, Radon finds that the quintet-triplet splitting change is almost linear as a function of the percentage of exact exchange used.While these results are qualitatively the same as the results presented here, the overall energy stabilization due to FLOSIC is only two thirds as large as that predicted by Radon.Since the results presented here exclude self-consistency, that was effectively noted to be important by Radon, we anticipate that future work which includes self-consistency, along the lines discussed by Yang et al. [53], may prove to be quantitatively important in predicting the precise triplet-quintet splitting.For a discussion of self-consistent treatments of transition-metal atoms within this formalism, see Kao et al. [54].

The S = 1 Porphyrin Structures
As discussed above, within GGA, there are two electronically distinct S = 1 spin configurations.In one case, the d xz /d yz minority states are half occupied with an unoccupied d 2 z minority orbital.In the other case, it is the other way around.The former state is Jahn-Teller unstable.These electronically distinct configurations may be found by starting potentials that favor ferromagnetic Fe-N alignment in one case and antiferromagnetic alignment in the other case.A starting potential that favors antiferromagnetic alignment affects the initial interactions between the occupied Fe(3d) and N(2p) states and allows for convergence to different electronic states.We have relaxed both of these geometries completely and determined that noticeable relaxation in the JT state is observed while the other maintains a D 4d symmetry.In Table 1, we simply describe the LUMO for the minimum-energy S = 1 state as "LUMO" since it does not have significant Fe(3d) character.All the other orbitals described in Table 1 showed localized Fe(3d) character and are labelled accordingly.For example, for the S = 1 (2q) state, the LUMO is mixed by d xz and d yz of iron, but the charges captured by the projected wavefunctions are only 0.29 e, which implies that the electrons are delocalized.

Zero-Point Vibrational Contributions and Vibrational Signatures
One possible way of distinguishing the low-spin and high-spin moments are to look for spin-dependent changes in the vibrational spectra that are strongly correlated with either Fe or N displacements.The softening the the high-spin zero-point energies can be observed in the infra-red (IR) and Raman modes.For example, near 190 cm −1 , the S = 2 configuration has two strong Raman modes, split by 10 cm −1 , and a pair of weak/moderate degenerate Raman modes that are hidden beneath the strong lower-energy mode (see lower panel of Figure 5).As shown in the upper panel of this figure, the splittings and shifts, of the strong modes, change significantly with the spin of the Fe(II).The weak/moderate modes shift upward by 15 cm −1 .The two strong modes shift upward by 30-60 cm −1 , respectively, depending on spin.Analysis of the vibrational modes shows that these modes have major contributions from the Fe and N atoms.The weak modes correspond to a pair of vibrational states associated with the nitrogen atoms vibrating out of plane.The lower-energy strong Raman mode is due to in-plane breathing of the N atoms with N atoms aligned along the y-axis breathing out of phase with the N atoms aligned along the x-axis.The slightly higher-energy mode is also due to the in-plane nitrogen atoms beating against the Fe but with a different phase.Calculation of the S = 1 Raman spectra, within PBE-GGA or PW92-LDA was complicated by the fact that there are two energetically similar S = 1 states.Application of electric fields seemed to mix these states, which made it difficult to reliably determine the derivatives of the polarizability tensors [31] needed to calculate the Raman spectra.Since the inclusion of the self-interaction correction significantly opens the HOMO/LUMO gap, it is reasonable to expect that this instability will not be present in a fully self-consistent SIC method.
As shown in Figure 6, the S = 2 geometry also has a pair of much weaker ungerade IR-active modes at 243 cm −1 .This pair of vibrational modes harden to 303 cm −1 in the S = 0 configuration.We have determined that these IR-active modes correspond to vibrations of the Fe ion along the (110) directions, to avoid the N atoms, in the plane of the porphyrin.These modes are expected to be sensitive to the presence of the high-energy d x 2 −y 2 states.As in the case of the gerade modes associated with N displacements, it is very reasonable to expect that this type of IR active mode could be significantly softened by a change in spin of the porphyrin.While these states are significantly changed by a change in the spin of the Fe ion, they are also very insensitive to changes in the density-functional or bending of the porphyrin.For example, the the PW92-LDA vibrational energies of the bent geometry agree with the PBE-GGA vibrational energies of the flat porphyrin to within 4 cm −1 .The vibrational zero-point energies presented in Table 2 correspond to the three states with HOMO/LUMO gaps.There is an additional (0.05 eV) stabilization of the S = 2 spin configuration relative to the S = 1 configuration when zero point energies are accounted for.This is due to the uniform softening observed in the IR and Raman spectra.

Conclusions
In this paper, we have applied a new SIC methodology, based upon Löwdin-orthonormalized Fermi Orbitals to the tetra-coordinated porphyrin molecule.Within the current implementation of the methodology, it is clear that the inclusion of this self-interaction-corrected formulation has significant effects on the energetics of undercoordinated Fe(II), as a function of spin.The results here suggest that the lowest energy spin state, within FLOSIC-PW92, is qualitatively different than the one predicted by PBE-GGA and PW92-LDA and that this difference appears to be in accordance with our analysis of experiments.It further appears that FLOSIC-PW92 prefers a non-planar S = 2 geometry, which also seems to be in accordance with experiments but in discordance with the vibrationally stable planar geometries predicted by the density-functional calculations performed here and by others.There are good reasons to expect that such corrections could be important in other transition-metal compounds if strong ligand-field effects are unable to control the d-state filling or in cases that ligands are missing.The results discussed here have not considered self-consistent relaxation of the KS spin density matrices.It is clear that this improvement should now be a high priority prior to additional tests of the method.However, overall, it would appear that the unitarily invariant inclusion of SIC provides relatively significant improvements of local density functionals.
In addition to using infra-red, Raman, and the high-field electron paramagnetic spectra, calculated here, comparison of calculated nuclear-magnetic resonance and Mössbauer spectra [59] would help significantly in determining whether the additional accuracy afforded by the self-interaction corrected methods will provide the chemical accuracy that is needed to fully understand the mechanistic details associated with metal-center catalysis, photocatalysis, spin-based quantum information, and gas-handling.

Figure 1 .
Figure1.Geometry of the S = 0 Fe(II)-porphyrin and its Fermi-orbital descriptors.The Fermi-orbital descriptors are presented in red.The central iron is in orange, the nitrogens are in blue, the carbons are in gray, and the hydrogens are in white.The hovering descriptors seem to be associated with a degree of aromaticity in this molecule.A discussion of the relationship between hovering FODs and aromaticity is provided in Refs.[32,42,49].

Figure 2 .
Figure 2. Orbital plot of the S = 2 majority-spin highest occupied molecular orbital, with d x 2 −y 2 symmetry in the bent geometry.This spin state and geometry is found to be the most stable by 0.055 eV.Within PW92-LDA, it is unstable by 1.24 eV.

Figure 3 .
Figure 3. Orbital plot of the S = 2 minority-spin highest occupied molecular orbital in the bent geometry.For this case, there is a slight Jahn-Teller distortion which preferentially occupies one of the otherwise degenerate d xz /d yz levels.The resulting d xz /d yz gap of the minority-spin down manifold is 0.28 eV.With inclusion of SIC, the gap between the occupied and unoccupied levels should increase significantly which makes the Jahn-Teller distortion even more favorable.

Figure 4 .
Figure 4.The orbital plot of the lowest energy S = 2 majority-spin d yz -orbital in the bent geometry.Upon bending, the d yz delocalizes slightly and takes on significant contributions from the pentagonal carbon cycles.With full inclusion of being self-consistent, this state is expected to localize in response to the coulomb SIC and further stabilize the S = 2 state relative to other spin states.

Figure 5 .Figure 6 .
Figure 5. Raman Spectra as a function of spin for the Fe(II)-porphyrin molecules.