New Insight into the Ground State of FePc: A Diffusion Monte Carlo Study

We have applied DMC to evaluate relative stability of the possible electronic configurations of an isolated FePc under D 4h symmetry, considering some fixed nodes generated from different methods. They predict A 2g ground state consistently, supporting preceding DFT studies, with confidence overcoming the ambiguity about exchange-correlation (XC) functionals. By comparing DMC with several XC, we clarified the importance of the short-range exchange to describe the relative stability. We examined why the predicted A 2g is excluded from possible ground states in the recent ligand field based model. Simplified assumptions made in the superposition model are identified to give unreasonably less energy gain for A 2g when compared with the reality. The state is found to have possible reasons for the stabilization, reducing the occupations from an unstable anti-bonding orbital, avoiding double occupation of a spatially localized orbital, and gaining exchange energy by putting a triplet spin pair in degenerate orbitals.


DFT calculations
We performed DFT calculations using Gaussian09 10 . We run all-electron calculations with def2QZVP basis set to compare DFT results. We used Burkatzki pseudo potentials with triple-ζ valence basis set to generate the trial nodes for DMC. 11 To get a symmetry-adapted state for each electronic configuration in a DFT simulation, we used "guess=alter" and "scf=symm" implemented in Gaussian09: We first give an initial guess appropriate for the target state by "guess=alter", and then fix the symmetries of the all occupied orbitals during its SCF procedure by "scf=symm".

DFT+U calculations
We have performed the DFT+U calculations for FePc molecule by using a simplified version of Cococcioni and de Gironcoli 12 , as implemented in QUANTUM ESPRESSO package 13 . Several different values (0, 2 and 4 eV) have been considered for Hubbard U parameter for Fe 3d orbitals. We have employed ultrasoft pseudopotentials generated with the Rappe-Rabe-Kaxiras-Joannopoulos recipe 14 to represent electron-ion interaction. The electronic exchange-correlation potential was calculated within the generalized gradient approximation (GGA) using the scheme of Perdew-Burke-Ernzerhof (PBE) 15 and the spin-polarization was taken into account. The electronic wave functions were expanded in plane waves with an energy cutoff of 35 Ry while for the charge density the energy cutoff was taken to 350 Ry. The isolated FePc molecule was simulated in a simple tetragonal cell of 27 × 27 × 12 Å 3 . Brillouin-zone integrations were approximated using a Γ point. The atomic positions of FePc molecule were optimized till the residual forces were less than 0.01 eV/Å.

Jastrow factor
We adopted a Jastrow factor 16 multiplied by determinant(s) to form a guiding function for DMC, imposing Kato's cusp conditions. 17 We used a function form implemented in CASINO 18 including electron-electron (u), electron-nuclei (χ), and electron-electron-nuclei ( f ) terms. Considering spin polarizations, u and χ ( f ) terms are expanded upto 8th (2nd) order of the power of inter-particle distances, getting total 144 variational parameters. Cutoff lengths for these terms are fixed as the recommended values by the implementation, and all the linear variational parameters are optimized by 'varmin-linjas' scheme. 19 For all-electron DMC, we also used the cusp-correction scheme 20 for possible electron-nuclei coalescence. All the present QMC simulations were done using CASINO (ver. 2.13) 18 .

Extrapolation of time step error
For time-step error corrections in DMC, 21 we used an extrapolation scheme 22 using two time steps, τ 1 = τ max and τ 2 = τ max /4. τ max is taken so that it is the largest possible below which the error is proportional to τ. From the results, E j ± σ j by τ j Figure S1. The extrapolation of DMC results. We can identify the ground state as A 2g with 1σ statistical confidence, while we cannot for other excited states. Except E g (a), linear extrapolations below τ max = 0.0004 works well. For E g (a), we extrapolate using τ max = 0.0002, instead.
For pseudo potential calculations, τ 1 = 4.0 × 10 −3 and τ 2 = 1.0 × 10 −3 were chosen for DFT-DMC. These values are larger than those of SA-CASCF-DMC. A necessary resolution in time step is larger in the pseudo potential case than in the all-electron because random walkers diffuse on shallower potentials.

Effects of geometry differences
As mentioned above, SA-CASSCF-DMC restricts us to use the same geometry 23 to all the states of electronic configurations. When the optimized geometry for each state largely differs from each other, this restriction could make the estimation poor. SA-CASSCF-DMC for an acrolein molecule 8 seems to be the case that it gave the overestimation of the excitation energy by ∼150 meV: Though it is not explicitly stated in their paper, 8 the bond length between carbon and oxygen gets elongated by 8 % 24, 25 when the system is excited. Because of the restriction, however, SA-CASSCF-DMC cannot take into account the relaxation energy gain by the elongation, and this is quite likely to be an origin of the overestimation. To examine if this matters in the present case, we evaluated the energy gains by the relaxation from the geometry used in SA-CASSCF, as tabulated in Tab. S1. It is confirmed that the gains remain within 1.3 %, corresponding to 7.7 meV which is negligibly small Table S1. The energy gains and the bond length owing to geometry optimization. B3LYP-DFT estimations of geometry relaxation gains, ∆E, and changes in the Fe-N bond lengths (R Fe−N ). ∆E is defined as the gain when the geometries are optimized from the common structure 23 used in the present SA-CASSCF. For each E g state, two bond lengths are given because it falls into D 2h from D 4h by the relaxation.
compared to the statistical errors and to the energy scale in Fig. 1 in the paper. The largest relaxation is found to occur on the bond between iron and neighboring nitrogen, which is confirmed to be within 0.26 % at most. The geometry insensitivity to the occupations to be considered is, incidentally, in accordance with a report 26 that the Fe-N bonding length is mainly dominated by the occupation number of d x 2 −y 2 , which is not considered here, though the conclusion is drawn from Fe porphyrin case. The insensitivity could justify the use of the same geometry to evaluate the relative stabilities among the states.