Hybrid Exchange-correlation Functionals Research Articles

Efficient real space formalism for hybrid density functionals.

We present an efficient real space formalism for hybrid exchange-correlation functionals in generalized Kohn-Sham density functional theory (DFT). In particular, we develop an efficient representation for any function of the real space finite-difference Laplacian matrix by leveraging its Kronecker product structure, thereby enabling the time to solution of associated linear systems to be highly competitive with the fast Fourier transform scheme while not imposing any restrictions on the boundary conditions. We implement this formalism for both the unscreened and range-separated variants of hybrid functionals. We verify its accuracy and efficiency through comparisons with established planewave codes for isolated as well as bulk systems. In particular, we demonstrate up to an order-of-magnitude speedup in time to solution for the real space method. We also apply the framework to study the structure of liquid water using abinitio molecular dynamics, where we find good agreement with the literature. Overall, the current formalism provides an avenue for efficient real-space DFT calculations with hybrid density functionals.

A simple method for calibration of hybrid exchange–correlation functionals for precise calculations of fluorescence wavelengths of dibenzoylmethanatoboron difluoride exciplexes with benzene, alkylbenzenes and pyridine

A simple approach to precise calculation of fluorescence wavelengths of charge-transfer exciplexes formed by dibenzoylmethanatoboron difluoride (DBMBF2) with benzene, alkylbenzenes and pyridine by varying the percentage of Hartree–Fock exchange energy (HFX) in the hybrid BHHLYP exchange–correlation functional is proposed. A correlation revealed between the optimal HFX parameter and the degree of charge transfer (CT) in the exciplexes provides a basis for the pre-assessment of the required HFX parameter in the BHHLYP functional from the calculated CT value in an exciplex of interest.

Topological nonsymmorphic insulator versus Dirac semimetal in KZnBi

KZnBi was discovered recently as a new three-dimensional Dirac semimetal with a pair of bulk Dirac fermions in contrast to the trivial insulator reported earlier. In order to address this discrepancy, we have performed electronic structure and topological state analysis of KZnBi using the local, semilocal, and hybrid exchange-correlation (XC) functionals within the density functional theory framework. We find that various XC functionals, including the SCAN meta-GGA and hybrid functional with 25% Hartree–Fock (HF) exchange (HSE06), resolve a topological nonsymmorphic insulator state with the glide-mirror protected hourglass surface Dirac fermions. By carefully tuning the XC strength in modified Becke-Johnson (mBJ) potential, we recover the correct orbital ordering and Dirac semimetal state of KZnBi. We further show that increasing the default HF exchange in hybrid functional (40%) can also capture the desired Dirac semimetal state with the correct orbital ordering of KZnBi. The calculated energy dispersion and carrier velocities of Dirac states are found to be in excellent agreement with the available experimental results. Our results demonstrate that KZnBi is a unique topological material where large XC effects are crucial to producing the Dirac semimetal state.

A quantum chemical investigation of the second hyperpolarizability of p-nitroaniline.

Recent measurements of the third harmonic scattering responses of molecules have given a new impetus for computing molecular second hyperpolarizabilities (γ) and for deducing structure-property relationships. This paper has employed a variety of wavefunction and density functional theory methods to evaluate the second hyperpolarizability of the p-nitroaniline prototypical push-pull π-conjugated molecule, addressing also numerical aspects, such as the selection of an integration grid and the impact of the order of differentiation vs the achievable accuracy by using the Romberg quadrature. The reliability of the different methods has been assessed by comparison to reference Coupled-Cluster Singles and Doubles with perturbative treatment of the Triples results. On the one hand, among wavefunction methods, the MP2 scheme offers the best accuracy/cost ratio for computing the static γ. On the other hand, using density functional theory, γ remains a challenging property to compute because all conventional, global hybrid or range-separated hybrid, exchange-correlation functionals underestimate static γ values by at least 15%. Even tuning the range-separating parameter to minimize the delocalization errors does not enable to improve the γ values. Nevertheless, the original double-hybrid B2-PLYP functional, which benefits from 27% of PT2 correlation and 53% Hartree-Fock exchange, provides accurate estimates of static γ values. Unfortunately, the best performing exchange-correlation functionals for γ are not necessarily reliable for the first hyperpolarizability, β, and vice versa. In fact, the β of p-nitroaniline (pNA) could be predicted, with a good accuracy, with several hybrid exchange-correlation functionals (including by tuning the range-separating parameter), but these systematically underestimate γ. As for γ, the MP2 wavefunction method remains the best compromise to evaluate the first hyperpolarizability of pNA at low computational cost.

Understanding the role of oxygen-vacancy defects in Cu2O(111) from first-principle calculations

The presence of defects, such as copper and oxygen vacancies, in cuprous oxide films determines their characteristic carrier conductivity and consequently their application as semiconducting systems. There are still open questions on the induced electronic re-distribution, including the formation of polarons. Indeed, to accurately reproduce the structural and electronic properties at the cuprous oxide surface, very large slab models and theoretical approaches that go beyond the standard generalized gradient corrected density functional theory are needed. In this work we investigate oxygen vacancies formed in proximity of a reconstructed Cu2O(111) surface, where the outermost unsaturated copper atoms are removed, thus forming non-stoichiometric surface layers with copper vacancies. We address simultaneously surface and bulk properties by modelling a thick and symmetric slab, to find that hybrid exchange-correlation functionals are needed to describe the oxygen vacancy in this system. Our simulations show that the formation of oxygen vacancies is favoured in the sub-surface layer. Moreover, the oxygen vacancy leads to a splitting and left-shift of the shallow hole states in the gap, which are associated with the deficiency of copper at the surface. These findings suggest that surface electronic structure and reactivity are sensitive to the presence of oxygen vacancies, also when the latter are formed deeper within the film.

Assessing the accuracy of hybrid exchange-correlation functionals for the density response of warm dense electrons.

We assess the accuracy of common hybrid exchange-correlation (XC) functionals (PBE0, PBE0-1/3, HSE06, HSE03, and B3LYP) within the Kohn-Sham density functional theory for the harmonically perturbed electron gas at parameters relevant for the challenging conditions of the warm dense matter. Generated by laser-induced compression and heating in the laboratory, the warm dense matter is a state of matter that also occurs in white dwarfs and planetary interiors. We consider both weak and strong degrees of density inhomogeneity induced by the external field at various wavenumbers. We perform an error analysis by comparing with the exact quantum Monte Carlo results. In the case of a weak perturbation, we report the static linear density response function and the static XC kernel at a metallic density for both the degenerate ground-state limit and for partial degeneracy at the electronic Fermi temperature. Overall, we observe an improvement in the density response when the PBE0, PBE0-1/3, HSE06, and HSE03 functionals are used, compared with the previously reported results for the PBE, PBEsol, local-density approximation, and AM05 functionals; B3LYP, on the other hand, does not perform well for the considered system. Additionally, the PBE0, PBE0-1/3, HSE06, and HSE03 functionals are more accurate for the density response properties than SCAN in the regime of partial degeneracy.

Dislocations as natural quantum wires in diamond

We study the electronic properties of the glide set of dislocations in diamond from first principles using hybrid exchange correlation functionals and find that the atomic-scale dislocation core states give rise to a prototypical one-dimensional (1D) band structure, i.e., natural quantum wires. The position and character of the core states varies strongly with local structure, where mixed dislocations with dangling bonds exhibit a 1D metallic band with a characteristic 1D density of states ($1/\sqrt{E}$). This 1D Fermi gas is spatially localized to a single atomic diameter orbital chain along the dislocation core. When the dangling bonds within the core are reconstructed, the 1D metallic band disappears. In contrast, pure edge dislocations in diamond reveal a 1D semiconductor with a direct band gap of 3.21 eV. These calculations provide a possible explanation to the long-standing observation of a blue luminescence band correlated with dislocations in diamond. This opens the door to using dislocations as 1D quantum phases with functional (electronic and optical) properties arising from the atomic-scale core states.

The structural, mechanical, thermal, electronic and optical properties of halide perovskites [formula omitted]: First-principles investigations

The structural, mechanical, thermal, electronic and optical properties of inorganic Pb-free double perovskites Cs2TiX6(X=Cl,Br,I) compounds are investigated. Using the first-principles with full-potential linearized augmented plane wave (FP-LAPW) method, supported the density functional theory, and therefore the generalized gradient approximation (GGA-PBE) and Hybrid exchange-correlation functionals with and without spin-orbit coupling. The computation of the lattice constants showed an honest agreement between our results and the experimental results and the theoretical data available. Using the three independent elastic constants C11, C12 and C44, various mechanical and thermal properties are obtained such as: bulk modules B, shear modules G, Young's modules Y, Pough's ratio B/G, Frantesvich ratio G/B, Poisson's ratio υ, Shear anisotropy character A and Debye temperature ΘD for Cs2TiI6, Cs2TiBr6 and Cs2TiCl6 respectively. Our calculated elastic constant results show that Cs2TiI6, Cs2TiBr6 and Cs2TiCl6 are mechanically stable, anisotropic and ductile in nature. Furthermore, the Cs2TiCl6 compound has the highest Debye temperature, which indicates the highest melting temperature. The electronic band structures reveal the semiconductor conductivity of the three compounds, with a band gap suitable for photovoltaic applications, within the range of 1.5–2.96 eV. Finally, Cs2TiCl6(X=Cl,Br,I) HPs compounds exhibit excellent optical properties like optical absorption, conductivity and refraction index.

Exciton dispersion and exciton-phonon interaction in solids by time-dependent density functional theory.

Understanding, predicting, and ultimately controlling exciton band structure and exciton dynamics are central to diverse chemical and materials problems. Here, we have developed a first-principles method to determine exciton dispersion and exciton-phonon interaction in semiconducting and insulating solids based on time-dependent density functional theory. The first-principles method is formulated in planewave bases and pseudopotentials and can be used to compute exciton band structures, exciton charge density, ionic forces, the non-adiabatic coupling matrix between excitonic states, and the exciton-phonon coupling matrix. Based on the spinor formulation, the method enables self-consistent noncollinear calculations to capture spin-orbital coupling. Hybrid exchange-correlation functionals are incorporated to deal with long-range electron-hole interactions in solids. A sub-Hilbert space approximation is introduced to reduce the computational cost without loss of accuracy. For validations, we have applied the method to compute the exciton band structure and exciton-phonon coupling strength in transition metal dichalcogenide monolayers; both agree very well with the previous GW-Bethe-Salpeter equation and experimental results. This development paves the way for accurate determinations of exciton dynamics in a wide range of solid-state materials.

Effects of Selenium Doping in Zinc Telluride from First Principles

Point defects and compensation in zinc telluride have been subjects of study for several decades due to the compound’s scientific/technological applications. Despite this, there remains some disagreement in the literature as to the fine details of its defect chemistry. Recently, it has also come under renewed interest due to the emerging Cd(1–x)ZnxTe(1–y)Sey material family, used in nuclear radiation detector applications. In this article, we use state-of-the-art computational tools to re-examine the native defect chemistry of ZnTe and explore how it changes upon addition of selenium using density functional theory with screened hybrid exchange–correlation functionals, thermochemical models of the equilibration environment, and statistical mechanics models of the point defect ensemble.

On native point defects in ZnSe

Aiming at a fundamental understanding of the defect chemistry of pure ZnSe for optical and quantum applications, systematic density functional theory calculations with hybrid exchange-correlation functionals were performed to build an accurate database of native defects in ZnSe, including isolated defects and first nearest-neighbor defect–defect complexes. From the defect formation energies, zinc vacancy is found to be the most prevalent defect as the Fermi level approaches the conduction band edge, while zinc interstitial in the selenium tetrahedron and selenium vacancy become the most prevalent defects as the Fermi level approaches the valence band maximum. The divacancy complex, consisting of first nearest-neighboring zinc and selenium vacancies, is also found to have a favorable binding energy across the entire bandgap. Its formation energy is, however, always higher than either the isolated zinc or selenium vacancy, meaning it will never be the predominant defect in equilibrium. Finally, a point defect with extended spin coherence in Fluorine-implanted ZnSe was recently discovered, and it was found to exhibit a broad emission peak centered at 2.28 eV. The identity of this defect was determined to be either zinc vacancy or its associated complex according to the electron paramagnetic resonance measurements. Explicit simulations of the optical signatures of all zinc vacancy-related native defects were conducted here, showing that both zinc vacancy and divacancy are the most likely native defect contributors to that peak.

KSSOLV 2.0: An efficient MATLAB toolbox for solving the Kohn-Sham equations with plane-wave basis set

KSSOLV (Kohn-Sham Solver) is a MATLAB toolbox for performing Kohn-Sham density functional theory (DFT) calculations with a plane-wave basis set. KSSOLV 2.0 preserves the design features of the original KSSOLV software to allow users and developers to easily set up a problem and perform ground-state calculations as well as to prototype and test new algorithms. Furthermore, it includes new functionalities such as new iterative diagonalization algorithms, k-point sampling for electron band structures, geometry optimization and advanced algorithms for performing DFT calculations with local, semi-local, and hybrid exchange-correlation functionals. It can be used to study the electronic structures of both molecules and solids. We describe these new capabilities in this work through a few use cases. We also demonstrate the numerical accuracy and computational efficiency of KSSOLV on a variety of examples. Program summaryProgram title: Kohn-Sham Solver 2.0 (KSSOLV 2.0)CPC Library link to program files:https://doi.org/10.17632/pp8vgvfcv4.1Developer's repository link:https://bitbucket.org/berkeleylab/kssolv2.0/src/release/Licensing provisions: BSD 3-clauseProgramming language:: MATLABNature of problem: KSSOLV2.0 is used to perform Kohn-Sham density functional theory based electronic structure calculations to study chemical and material properties of molecules and solids. The key problem to be solved is a constrained energy minimization problem, which can also be formulated as a nonlinear eigenvalue problem.Solution method: The KSSOLV 2.0 implements both the self-consistent field (SCF) iteration with a variety of acceleration strategies and a direct constrained minimization algorithms. It is written completely in MATLAB and uses MATLAB's object oriented programming features to make it easy to use and modify.

DensityTool: A post-processing tool for space- and spin-resolved density of states from VASP

The knowledge of the local electronic structure of heterogeneous solid materials is crucial for understanding their electronic, magnetic, transport, optical, and other properties. VASP, one of the mostly used packages for density-functional calculations, provides local electronic structure either by projecting the electronic wave functions on atomic spheres, or as a band-decomposed partial charge density. Here, we present a simple tool which takes the partial charge density and the energy eigenvalues calculated by VASP as input and constructs local charge and spin densities. The new data provides a much better spatial understanding than the projection on the atomic spheres. It can be visualized directly in the real space e.g. with Vesta, or averaged along planes spanned by two of the lattice vectors of the periodic unit cell. The plane-averaged local (spin) density of states can be easily plotted e.g. as color-coded data using almost any plotting program. DensityTool can be applied to visualize and understand the local electronic structure of any system calculated with VASP. We expect it to be useful especially for researchers concerned with inhomogeneous systems, such as interfaces, defects, surfaces, adsorbed molecules, or hybrid inorganic-organic composites. Program summaryProgram Title: DensityToolCPC Library link to program files:https://doi.org/10.17632/26mcg4jwvd.1Developer's repository link:https://github.com/llodeiro/DensityToolLicensing provisions: MITProgramming language: FORTRANSupplementary material: A complementary User Manual and examples can be found in the developer's repository link.Nature of problem: Total and local density of states (LDOS) are widely-used quantities for understanding the electronic structure of condensed matter. In VASP, one of the mostly used packages for numerically solving the Kohn-Sham equations of density-functional theory, one does not directly obtain the LDOS. Usually, one works with partial density of states (PDOS), where the wave function is projected onto atomic orbitals. Alternatively, it is possible to calculate the band- and wavevector-dependent charge densities. Together with the band structure, LDOS can be obtained using post-processing tools.Solution method: DensityTool combines the band structure with band- and wavevector-dependent charge densities calculated by VASP to construct LDOS. In addition, spin-resolved LDOS (LSDOS) can be calculated for magnetic systems. Thus obtained LDOS and LSDOS can be further processed with DensityTool. In particular, it can be averaged along chosen spatial directions, which is convenient for example to demonstrate the spatially-resolved electronic structure of inhomogeneous systems.DensityTool can handle all three types of approaches implemented in VASP (non-magnetic, collinear magnetic, non-collinear magnetic including spin-orbit coupling), as well as all types of self-consistent band calculations as DFT (e.g., LDA, GGA, mGGA, BJ, and hybrid exchange-correlation functionals) and HF.For convenience, the output of DensityTool follows the format of VASP (the CHGCAR file), which can be easily used in other programs for visualization. We provide a sample python script to plot averaged L(S)DOS from the calculated data.

Length-Gauge Optical Matrix Elements in WIEN2k

Hybrid exchange-correlation functionals provide superior electronic structure and optical properties of semiconductors or insulators as compared to semilocal exchange-correlation potentials due to admixing a portion of the non-local exact exchange potential from a Hartree–Fock theory. Since the non-local potential does not commute with the position operator, the momentum matrix elements do not fully capture the oscillator strength, while the length-gauge velocity matrix elements do. So far, length-gauge velocity matrix elements were not accessible in the all-electron full-potential WIEN2k package. We demonstrate the feasibility of computing length-gauge matrix elements in WIEN2k for a hybrid exchange-correlation functional based on a finite difference approach. To illustrate the implementation we determined matrix elements for optical transitions between the conduction and valence bands in GaAs, GaN, (CH3NH3)PbI3 and a monolayer MoS2. The non-locality of the Hartree–Fock exact exchange potential leads to a strong enhancement of the oscillator strength as noticed recently in calculations employing pseudopotentials (Laurien and Rubel: arXiv:2111.14772 (2021)). We obtained an analytical expression for the enhancement factor for the difference in eigenvalues not captured by the kinetic energy. It is expected that these results can also be extended to other non-local potentials, e.g., a many-body GW approximation.

The role of tagging atoms on the thermal stability and vibrational behavior of Nb9 clusters

Extensive comparative study of density functional theory calculations dedicated to analyzed the influence of the so-called tagging atoms on the thermal stability and vibrational behavior of small Nb9—Arn (n = 1–4) clusters are presented. We consider two-hybrid functional (HSEH1PBE and CAM-B3LYP) models to obtain minimum energy structures and, the relative concentrations of several Ar complexes configurations when varying the temperature of the system are determined by constructing partition functions and enthalpies. Under these considerations, it was found that the coexistence of two different isomers of the Nb9 clusters in a sample is highly probable. Additionally, some interesting points can be highlighted, for instance, i) the different hybrid exchange-correlation functionals modify the proportions of mole fractions in all cases, ii) the weighted average of the IR spectra proposed by the mole fractions have better agreement with the experimental data, iii) tagging atoms can induce significant variations in the infrared spectra, such variations consist mainly of noticeable changes in the relative intensity of the active infrared frequencies. In conclusion, even if the percentage of isomers shown in the mole fractions changes as a function of the hybrid exchange-correlation functional, the influence of tagging atoms makes some changes in the calculation of the partition function and in the simulated IR spectra. For this reason, it is suggested to include tagging atoms in theoretical calculations to have a better understanding of the experimental data.

Massively parallel linear-scaling Hartree-Fock exchange and hybrid exchange-correlation functionals with plane wave basis set accuracy.

We extend our linear-scaling approach for the calculation of Hartree-Fock exchange energy using localized in situ optimized orbitals [Dziedzic et al., J. Chem. Phys. 139, 214103 (2013)] to leverage massive parallelism. Our approach has been implemented in the onetep (Order-N Electronic Total Energy Package) density functional theory framework, which employs a basis of non-orthogonal generalized Wannier functions (NGWFs) to achieve linear scaling with system size while retaining controllable near-complete-basis-set accuracy. For the calculation of Hartree-Fock exchange, we use a resolution-of-identity approach, where an auxiliary basis set of truncated spherical waves is used to fit products of NGWFs. The fact that the electrostatic potential of spherical waves (SWs) is known analytically, combined with the use of a distance-based cutoff for exchange interactions, leads to a calculation cost that scales linearly with the system size. Our new implementation, which we describe in detail, combines distributed memory parallelism (using the message passing interface) with shared memory parallelism (OpenMP threads) to efficiently utilize numbers of central processing unit cores comparable to, or exceeding, the number of atoms in the system. We show how the use of multiple time-memory trade-offs substantially increases performance, enabling our approach to achieve superlinear strong parallel scaling in many cases and excellent, although sublinear, parallel scaling otherwise. We demonstrate that in scenarios with low available memory, which preclude or limit the use of time-memory trade-offs, the performance degradation of our algorithm is graceful. We show that, crucially, linear scaling with system size is maintained in all cases. We demonstrate the practicability of our approach by performing a set of fully converged production calculations with a hybrid functional on large imogolite nanotubes up to over 1400 atoms. We finish with a brief study of how the employed approximations (exchange cutoff and the quality of the SW basis) affect the calculation walltime and the accuracy of the obtained results.

Native oxide reconstructions on AlN and GaN (0001) surfaces

Properties of AlN/GaN surfaces are important for realizing the tunability of devices, as the presence of surface states contributes to Fermi level pinning. This pinning can influence the performance of high electron mobility transistors and is also important for passivation of the surface when developing high-power electronic devices. It is widely understood that both AlN and GaN surfaces oxidize. Since there are many possible reconstructions for each surface, it is a challenge to identify the relevant surface reconstructions in advance of a detailed simulation. Because of this, different approaches are often employed to down select initial structures to reduce the computational load. These approaches usually rely on either electron counting rules or oxide stoichiometry, as both of these models tend to lead to structures that are energetically favorable. Here we explore models from these approaches but also explore a reconstruction of the (0001) surface directly observed using scanning transmission electron microscopy with predictive density functional theory simulations. Two compositions of the observed surface reconstruction—one which obeys oxide stoichiometry and one which is cation deficient and obeys electron counting—are compared to reconstructions from the previous work. Furthermore, surface states are directly calculated using hybrid exchange-correlation functionals that correct for the underestimation of the bandgaps in AlN and GaN and improve the predicted positions of surface states within the gap. It is found that cation deficiency in the observed reconstruction yields surface states consistent with the experiment. Based on all of these results, we provide insight into the observed properties of oxidized AlGaN surfaces.

Hybrid functional analysis of planar metal-organic frameworks

The electronic energy band structure and optical properties of the planar metal-organic frameworks, namely M-HAB and M3[HITP] (M = Ni, Cu; HAB = hexaaminobenzene; HITP = hexaiminotriphenylene) were investigated in this work using two hybrid exchange-correlation functionals. The materials studied here contain nickel and copper, 3d electrons of which should be considered as strongly correlated subsystem. These semi-core electrons are treated here by means of the hybrid HSE06 and PBE0 functionals. We have found that Ni-HAB and Cu-HAB are metals, whereas Ni3[HITP]2 and Cu3[HITP]2 are degenerate p-type semiconductors.

Ground-to-excited derivative couplings for the density functional-based tight-binding method: semi-local and long-range corrected formulations

A derivation of non-adiabatic coupling vectors for the density functional-based tight-binding method (DFTB) between ground and excited states is presented. The analytical result is valid both for semi-local and long-range corrected DFTB and includes all required Pulay terms. Electron-translation factors lead to a conceptual simplification of the Slater-Koster scheme for precomputed integrals. Compared to scalar couplings obtained from numerical derivatives, the present approach is computationally more efficient and can be applied to systems with hundreds of atoms. The accuracy of DFTB derivative couplings is assessed by comparison to full density functional theory (DFT) calculations using semi-local and hybrid exchange-correlation functionals with promising results. As exemplified by a case study of furan, DFTB provides non-adiabatic coupling vectors that are close to DFT counterparts in size and direction also in the vicinity of conical intersections.

Accurate Prediction of Mössbauer Hyperfine Parameters in Bis-Axially Coordinated Iron(II) Phthalocyanines Using Density Functional Theory Calculations: A Story of a Single Orbital Revealed by Natural Bond Orbital Analysis.

Density Functional Theory (DFT) calculations coupled with several exchange-correlation functionals were used for the prediction of Mössbauer hyperfine parameters of 36 bis-axially coordinated iron(II) phthalocyanine complexes with the general formulas PcFeL2, PcFeL'L″, and [PcFeX2]2-, including four new compounds. Both gas-phase and PCM calculations using BPW91 and MN12L exchange-correlation functionals were found to accurately predict both Mössbauer quadrupole splittings and the correct trends in experimentally observed isomer shifts. In comparison, hybrid exchange-correlation functionals underestimated quadrupole splittings, while still accurately predicted isomer shifts. Out of ∼40 exchange-correlation functionals tested, only MN12L was found to correctly reproduce quadrupole splitting trends in the PcFeL2 complexes coordinated with phosphorus-donor axial ligands (i.e., P(OnBu)3 ≈ P(OEt)3 < PMe3 < P[(CH2O)2CH2]-p-C6H4NO2 < PEt3 ≈ PnBu3). Natural Bond Orbital (NBO) analysis was successfully used to explain the general trends in the observed quadrupole splitting for all compounds of interest. In particular, the general trends in the quadrupole splitting correlate well with the axial ligand dependent, NBO-predicted population of the 3dz2 orbital of the Fe ion and are reflective of the hypothesis proposed by Ohya and co-workers ( Inorg. Chem., 1984, 23, 1303) on the adaptability of the phthalocyanine's π-system toward Fe-Lax interactions. The first X-ray crystal structure of a PcFeL2 complex with axial phosphine ligands is also reported.