All the ground and excited state energies and wave function calculations are performed on the experimentally reported X-ray crystal structures' geometries. The molecular coordinate system is chosen in such a way that Fe(I) is the origin of the coordinate system and the z-axis points approximately towards the ligands. The correlated calculations are done using Complete Active Space Self Consisent Field (CASSCF) together with N-Electron Valence Perturbation Theory (NEVPT2). The Fe(I) possess d 7 electronic configuration giving rise to 10 quartets (S=3/2) and 40 doublets (S=1/2) states. To benchmark the methodology, we performed the calculation for complex 1 considering all the roots of quartet and doublet and another calculation with all roots of quartet only. However, the doublet states do not show any significant contribution to the D values, therefore, the calculations are performed accounting for only 10 quartet states (see Table in SI). The CASSCF energy levels are obtained by state-averaging these states in the active space consisting of 7 electrons distributed in 5 3d orbitals i.e., CAS (7,5). The effect of dynamical correlations is included by performing NEVPT2 calculations on top of CASSCF converged wavefunction. Scalar relativistic effects are taken into account employing second-order Douglas-Kroll-Hess triple-DKH-def2-TZVP basis set and auxiliary basis set is generated automatically using AUTOAUX keyword. Tight self-consistent field (SCF) convergence criteria (10 -10 au) is used in all the calculations. The spin-orbit interactions along with the zeeman interactions are incorporated by the quasi-degenerate perturbation theory (QDPT) approach. Magnetic parameters (D, E/D and g-tensors) have been computed using effective Hamiltonian theory as implemented in the ORCA (4.0.1.2) program package. Ab initio ligand field theory (AILFT) analysis based on wavefunction obtained from SA-CASSCF/NEVPT2 has been employed to obtain precise d-orbital splitting pattern of the studied complexes. Ab initio blocking barriers for relaxation mechanisms have been computed from SINGLE_ANISO module as implemented in Orca v5.0.1. The ab initio CASSCF and CASPT2 calculations are also performed in Open-Molcas (19.11). For CASSCF calculations in OpenMolcas, the relativistic contracted atomic natural orbital type basis sets i.e., [ANO-RCC...5s4p2d] for Fe, [ANO-RCC...4s3p] for Si, [ANO-RCC...3s2p] for N and C and [ANO-RCC...2s] for H are employed. To account for spin-orbit effects, AMFI spin-orbit operator is used and DKH Hamiltonian is considered for scalar relativistic effect. The AMFI (atomic mean field integral) spin-orbit operator is introduced to account for the spin-orbit effects. The scalar relativistic effect is considered using the DKH Hamiltonian. The resolution-of-identity approximation is employed to accelerate the calculation.