Influence of Electron Correlation on the Electronic Structure and

Article pubs.acs.org/JCTC

Influence of Electron Correlation on the Electronic Structure and Magnetism of Transition-Metal Phthalocyanines Iulia Emilia Brumboiu, Soumyajyoti Haldar, Johann Lüder, Olle Eriksson, Heike C. Herper, Barbara Brena,* and Biplab Sanyal Department of Physics and Astronomy, Uppsala University, SE-75120 Uppsala, Sweden S Supporting Information *

ABSTRACT: There exists an extensive literature on the electronic structure of transition-metal phthalocyanines (TMPcs), either as single molecules or adsorbed on surfaces, where explicit intra-atomic Coulomb interactions of the strongly correlated orbitals are included in the form of a Hubbard U term. The choice of U is, to a large extent, based solely on previous values reported in the literature for similar systems. Here, we provide a systematic analysis of the influence of electron correlation on the electronic structure and magnetism of several TMPcs (MnPc, FePc, CoPc, NiPc, and CuPc). By comparing calculated results to valence-band photoelectron spectroscopy measurements, and by determining the Hubbard term from linear response, we show that the choice of U is not as straightforward and can be different for each different TMPc. This, in turn, highlights the importance of individually estimating the value of U for each system before performing any further analysis and shows how this value can influence the final results. photovoltaics3 to molecular electronics and spintronics.4−8 It is important to note that the majority of these applications involve the adsorption of the TMPc molecules on metal surfaces or between metal electrodes. In addition, the interesting magnetic (spin state) and electronic properties (band gap, hybridization with ligand states, highest occupied molecular orbital (HOMO) character) of TMPcs originate essentially in the d-levels of the central metal atom and are influenced both by the ligands and by the presence of the surface/electrodes.9−14 Thus, the correct description of the dstates both in the presence and in the absence of a surface is crucial. Density functional theory (DFT) within the local density approximation (LDA) or the generalized gradient approximation (GGA), although successful and computationally efficient in describing basic structural and chemical properties of molecules and extended systems,15,16 is not sufficient for describing localized and strongly correlated electrons.17 Hybrid functionals, which involve the addition of a portion of exact exchange to the exchange and correlation functional, instead provide rather reliable results in what concerns the electronic structure of single molecules containing transition metals. Many studies on the valence-band (VB) electronic structure of TMPcs and transition-metal porphyrins show that hybrids provide a good agreement with the VB photoelectron spectra (PES).12,18−23 Unfortunately, the use of hybrid functionals for extended systems is computationally very expensive, because of the necessity of calculating the Hartree−Fock exact ex-

1. INTRODUCTION Metal phthalocyanines are highly symmetric organic molecules that incorporate a metal atom at the center of the organic ligand (as shown in Figure 1). Used as dyes on a regular basis,1,2 transition-metal phthalocyanines (TMPcs) have been the focus of scientific research for a long time, because of application prospects in a wide number of fields, ranging from

Figure 1. Molecular structure of a transition-metal phthalocyanine (TMPc). The two types of nitrogen atoms of the organic ligand are marked as Niso (directly bound to the metal) and Naza (N atoms that form the aza-bridge). The x- and y-axes are marked with arrows. The z-axis is perpendicular to the molecular plane. © XXXX American Chemical Society

Received: February 28, 2016

A

DOI: 10.1021/acs.jctc.6b00091 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation change.24,25 Moreover, it has been shown that the hybrids fail to describe correctly the electronic structure and magnetic properties of metallic surfaces.26 A solution to the problem of describing transition-metal molecules adsorbed on surfaces is the use of GGA augmented by a Hubbard U term. In this approach, an additional term is selectively added to the density functional theory (DFT) energy functional (eq 1), directly affecting only the localized and correlated d-orbitals of the transition metal in the system.27−29 The total energy is a functional of the charge density (ρσ) and the occupation numbers ({nIσ}) and can be written as follows:

values between 1 and 5 eV also have been reported.11,21,39,40 With regard to the other three phthalocyanines (CoPc, NiPc, and CuPc), the value Ueff = 3 eV is widely used,32,33,37,38,40,41 but a value of Ueff = 5 eV also has been employed.33,42,43 In the case of CuPc, a value of Ueff = 8.9 eV calculated from linear response has been reported by Mattioli and co-workers.9 Generally, the choice of U is based on previous reported values for similar systems, but a thorough analysis of the influence of the U parameter on the electronic and structural properties of these materials has not yet been performed. By comparing DFT +U calculations for various values of Ueff with hybrid functional results and PES measurements, our study aims to provide insight into how the molecular properties vary with Ueff and discuss which value could be appropriate for each single molecule. In addition, we compute the Hubbard U value, using a linear response approach29 (note that J is not calculated). Since the U term is defined as the second derivative of the energy, with respect to on-site occupation numbers,44 its value can be computed by studying the response of the system to a perturbation on the localized d-states.29,31

E DFT + U [ρ σ (r), {n Iσ }] = E DFT[ρ σ (r)] + E Hub[{n Iσ }] − Edc[{n Iσ }]

(1)

where EDFT is the DFT energy functional, EHub the Hubbard functional, defined for the localized orbitals, and Edc a “doublecounting” term that accounts for the contribution of the delectrons already included in EDFT. EDFT is a functional of the charge density, while EHub and Edc are functionals of the on-site occupation numbers. The Hubbard functional and the “double-counting” energy are expressed in terms of the effective Coulomb (U) and exchange (J) interactions, which, for d-electrons, are U=

J=

1 (2l + 1)2

U=

U = χ0−1 − χ −1

∑ ⟨m , m′|Vee|m , m′⟩ (2)

m,m′

1 2l(2l + 1)

∑

(3)

Here, Vee represents the electron−electron interaction and the states designated by the index m are the localized basis functions. There are two main approaches for implementing the DFT +U method practically. In the formulation proposed by Lichtenstein et al.,27 the Hubbard functional is rotationally invariant and contains fully orbital-dependent electronic interactions.27,29 The functionals EHub and Edc are defined as described in detail in refs 27, 29, 30, and 31. In a simplified approach, proposed by Dudarev et al.,28,29 the orbital-dependent electronic interactions are replaced by a spherical average. In this way, the rotational invariance is maintained, while the off-diagonal elements of the occupation matrix can be disregarded, making the computation of the Hubbard functional simpler. The exchange correction (J) is directly included in the Coulomb parameter and the effective value is defined as:

Ueff = U − J

(4)

⎛ Ueff ⎞ Iσ Iσ ⎟Tr[n (1 − n )] 2 ⎠

∑ ⎜⎝ I ,σ

(7)

2. COMPUTATIONAL METHODS We analyzed the electronic, magnetic, and structural properties of five TMPcs: MnPc, FePc, CoPc, NiPc, and CuPc. Several types of DFT calculations have been performed: hybrid functional calculations, using the hybrids B3LYP45 and HSE06;46−52 pure GGA calculations using the Perdew, Burke, and Ernzerhof (PBE) functional;53,54 GGA+U in the Dudarev implementation,28,29 using five different values of Ueff (from 1 eV to 5 eV) for each molecule; and a linear response calculation of U (eqs 6 and 7).29 In addition, in the case of CoPc, the vibrational frequencies have been calculated at two values of Ueff (2 and 4 eV) and compared to B3LYP results obtained using the Gaussian 09 D.01 quantum chemistry software (G09).55 The hybrid functional spin-polarized calculations were performed using G09. The hybrid functional by Becke45 and the screened hybrid functional HSE0646−52 were used in combination with the 6-31G(d,p) basis set56 for the C, N, and H atoms and the triple-ζ cc-pVTZ basis set57 for the metal atom. The metal d projected density of states was determined using the c2 method.58,59 The bar graphs corresponding to the discrete d-levels were broadened using Gaussian functions with a constant full width at half-maximum (fwhm) of 0.2 eV. The vibrational frequencies were calculated for CoPc using the geometries optimized at the same level of theory (i.e., B3LYP) and tight convergence criteria. The infrared (IR) intensities and Raman activities were determined as well. The Raman intensities (I(νi)) were calculated from the Raman activities (Si) using the following equation:60

This relationship is used in the expression of the term to be added to EDFT: E Hub − Edc =

(6)

where χ0 is the bare response matrix and χ is the self-consistent response matrix. The matrix elements are defined by the variation of the occupation number at site I (nI), with respect to an applied perturbation at site I′ (αI′), determined at the first iteration, in the case of χ0, and after self-consistency has been achieved, in the case of χ.31

⟨m , m′|Vee|m′, m⟩

m≠m′,m′

d2EDFT d(n I )2

(5)

Iσ

where n represent the occupation matrices for the atomic site I and spin σ. Many studies of TMPcs adsorbed on surfaces make use of DFT in its GGA implementation. Among these, several adopt the Hubbard U approach, most often in the formulation by Dudarev and co-workers.28 For MnPc, the Ueff value reported in the literature is 3 eV.32−34 Similarly, in the case of FePc, the majority of studies use a value of Ueff = 3 eV.21,33−38 Other B

DOI: 10.1021/acs.jctc.6b00091 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 2. (a) HOMO−LUMO gaps calculated using hybrid functionals represented as a function of the atomic number of the central metal (the symbols represented with full lines are B3LYP results, while those represented with dotted lines are determined by HSE06). (b) The variation in calculated HOMO−LUMO gaps with Ueff. (c) The relative change in the magnetic moments per metal atom as a function of Ueff.

I(νi) =

f ·Si(ν0 − νi)4

νi⎡⎣1 − e−hcνi /(kBT )⎤⎦

and correlation, in combination with ultrasoft pseudopotentials (from http://www.quantum-espresso.org) for the different atomic species.68 The cell size was the same as in the case of the VASP calculations, while the energy cutoff for the wave functions was 30 Ry and the kinetic energy cutoff for the charge density and potential was 360 Ry.

(8)

where I(νi) is the Raman intensity at the particular vibrational frequency νi (expressed in cm−1), f is a scaling factor (in our case, 1.0), Si are the computed Raman activities, and ν0 is the wavenumber corresponding to the laser wavelength used to excite the vibrational modes (632.8 nm, see Figure 3 and ref 81); h is Planck’s constant, c is the speed of light, and kB is the Boltzmann constant. T represents the temperature (293.15 K). The calculated IR and Raman spectra were normalized, with respect to the maximum intensity peak. A Gaussian broadening with fwhm = 10 cm−1 was added to all the calculated bar graphs in order to facilitate the comparison to experimental spectra. A similar broadening was added to the bar graphs corresponding to the normal modes, such that a comparison to the DFT+U results can be performed. No corrections were added to the calculated frequency scale. The GGA and GGA+U spin-polarized calculations were performed using the PBE53,54 exchange and correlation functional, as implemented in the Vienna ab initio simulation package (VASP).61−64 The projector augmented-wave (PAW) potentials65,66 were used to perform Γ-point calculations using a large cell (25 Å × 25 Å × 25 Å), with a minimum of 12 Å of vacuum separating the molecules in neighboring unit cells and an energy cutoff of 400 eV. The GGA+U calculations were performed in the Dudarev approach,28 using a fixed value for the J-parameter of 1.0 eV and varying the U parameter to obtain effective values of 1, 2, 3, 4, and 5 eV. All the systems were relaxed without symmetry constraints and the optimized structures were then used to statically determine the total and projected densities of states. In the case of CoPc, the vibrational frequencies at two different values of Ueff (2 and 4 eV) were determined using the finite differences approach without symmetry constraints. For these calculations, the geometries were first relaxed at the particular value of Ueff imposing tight convergence criteria for the forces (0.002 eV/Å). A Gaussian broadening of 10 cm−1 fwhm was added to the vibrational frequency bar graphs in order to facilitate the comparison between the DFT+U results and the G09 calculations. Finally, a linear response calculation as proposed by Cococcioni and De Gironcoli29,31 was performed on the VASP optimized structures, in order to determine the value of U. This calculation was carried out with the Quantum Espresso ab initio software,67 using the PBE functional53,54 for exchange

3. RESULTS AND DISCUSSION In the following, the influence of Ueff on the magnetic, electronic, and structural properties of the five TMPcs is thoroughly analyzed. First, we describe the variations of the magnetic moment per metal atom, of the HOMO−LUMO gaps, and of the metal−nitrogen bond lengths, followed by a paragraph on the vibrational frequencies calculated for CoPc. We then analyze the 3d electronic configurations of each phthalocyanine molecule for each value of Ueff. Comparisons with hybrid calculations and PES experiments are finally presented in the last part of this section. 3.1. HOMO−LUMO Gaps and Magnetic Moments. Figure 2 summarizes the variation in HOMO−LUMO gaps (Figures 2a and 2b) and calculated magnetic moments per metal atom (Figure 2c) with the change of Ueff. The magnetic moment is expressed as a percentage relative to the moment, calculated using a value of Ueff = 0 eV (Δμ (%) = [(μU − μ0)/ μ0] × 100). As can be seen in Figure 2b, the HOMO−LUMO gaps increase as Ueff increases, in the case of FePc, CoPc, and CuPc; however, they remain essentially unchanged in the case of MnPc and NiPc. It is interesting to note that the increase in the gaps of FePc and CoPc is initially strong, up to Ueff = 3 eV. Above this value, the curves flatten out and the increment of the HOMO−LUMO gaps becomes very small. The reason for this is the fact that, at low Ueff values, the character of the HOMO and LUMO orbitals of FePc and the character of the LUMO in CoPc are provided by the transition-metal atom. As the Ueff term is increased, the metal-contributed occupied (unoccupied) states move toward higher (lower) binding energies. Above the Ueff = 3 eV threshold, the frontier orbitals are defined instead by ligand states (C 2p, HOMO; C 2p and N 2p, LUMO), which are unaffected by U, and, therefore, the increase in the values of the HOMO−LUMO gaps stops. In the case of NiPc, the HOMO−LUMO gap is defined from the beginning by ligand states and thus remains essentially unchanged as Ueff is increased. Interestingly, in the case of MnPc, the HOMO is a hybrid molecular orbital containing C

DOI: 10.1021/acs.jctc.6b00091 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Table 1. Calculated M−Niso Bond Lengths Using the B3LYP and HSE06 Hybrid Functionals, in Comparison to the Different Ueff Resultsa

a

Mn−Niso (Å)

Fe−Niso (Å)

Co−Niso (Å)

Ni−Niso (Å)

Cu−Niso (Å)

exp

1.93 (ref 71)

1.93 (ref 72)

1.91 (ref 71)

1.87−1.91 (refs 73, 74)

1.93−1.95 (refs 73, 75−77)

B3LYP

1.96 1.95

1.95

1.94

1.92

1.97

HSE06

1.95 1.94

1.94

1.92

1.91

1.95

Ueff = 0 eV

1.95 1.94

1.94

1.92

1.91

1.96

Ueff = 1 eV

1.96 1.95

1.94

1.92

1.91

1.96

Ueff = 2 eV

1.97 1.95

1.94

1.93

1.91

1.96

Ueff = 3 eV

1.97 1.95

1.95

1.93

1.91

1.96

Ueff = 4 eV

1.97 1.95

1.95

1.93

1.91

1.96

Ueff = 5 eV

1.97 1.95

1.96

1.94

1.92

1.96

Functionals that give one value indicate that the D4h symmetry is preserved, while functionals that give two values indicate D2h symmetry.

spin-up and spin-down metal d-levels (see Figure S3 in the Supporting Information). More specifically, a small percentage of the spin-down unoccupied levels migrate to the occupied spin-up channel. The most affected d-orbital is the dxz (contributing to the HOMO), for which the spin-down occupation changes from ∼0.4 to ∼0.3 as Ueff increases to 5 eV. The increased moment on the Mn atom is compensated by an opposite moment on the ligand and the total magnetization of the molecule remains constant. The atoms that are the most involved are the isoindole nitrogens (Niso) and the alpha carbons (Cα), which both contribute with comparable amounts to the ligand moment. In the case of CuPc, the magnetic moment on Cu changes from 0.55 μB at Ueff = 0 eV to 0.62 μB at Ueff = 5 eV. The total magnetization of 1 μB is achieved as a result of components localized on the ligand isoindole N atoms (Niso) that sum to 0.35−0.3 μB as Ueff is increased from 0 to 5 eV. The remaining moment is due to small contributions from the rest of the ligand. In the hybrid functional calculations, Mulliken spin densities have been determined for each metal atom as the difference between the spin-up and spin-down Mulliken charges.69 The unpaired spins corresponding to the Mulliken spin densities can be assigned to magnetic moments and compared to those calculated by DFT+U. The hybrid functional moments for MnPc and CuPc approach values determined using the Hubbard U. Specifically, the B3LYP result of 3.43 μB for MnPc is very similar to the one calculated with Ueff = 3 eV (3.42 μB). The HSE06 value (3.57 μB) comes, instead, closer to the 3.54 μB calculated using Ueff = 4 eV. In the case of CuPc, the magnetic moments obtained by Ueff = 2 eV (0.57 μB) and

both metal (dπ) and C 2p contributions in equal amounts (B3LYP result22). The metal dπ-states remain pinned, despite the increase in Ueff. The other occupied d-orbitals, which do not contribute to the HOMO, move toward higher binding energies, as was the case for the other TMPcs. The pinning of the dπ is the reason for the lack of variation in the HOMO− LUMO gap with Ueff for MnPc. Finally, for CuPc, the variation in the gap has a low slope and appears linear. In CuPc, at low values of Ueff, the HOMO is characterized by both ligand and metallic contributions. On the other hand, the LUMO is entirely constituted by the metal dx2−y2 orbital. As Ueff is increased, the ligand HOMO remains in place, while the metal LUMO moves toward lower binding energies, giving rise to the observed trend in Figure 2b. We will further analyze these behaviors in the following part of this section, when the changes with Ueff of the projected densities of states (PDOS) will be described, in comparison to the total DOS. The GGA+U HOMO−LUMO gaps are consistently lower, compared to the hybrid functionals results (Figure 2a), by as much as 1.5 eV (MnPc, B3LYP). HSE06 and B3LYP are also different from each other, as B3LYP gaps are ∼0.5 eV larger in all five TMPcs. This difference is likely due to the different percentage of exact exchange in the two functionals. The changes in magnetic moments per metal atom, as a function of Ueff, are represented in Figure 2c. Variations are most notable in the cases of CuPc and MnPc, reaching 12% in the first case and 20% in the second, as Ueff increases to 5 eV. In contrast, the moment steadily, but only slightly increases (up to 3%) in the case of FePc and CoPc. NiPc (not shown in the figure) is predicted as nonmagnetic by all GGA+U calculations. The large variation in the magnetic moment on the Mn atom in MnPc is explained by the change in the occupation of the D

DOI: 10.1021/acs.jctc.6b00091 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

Figure 3. IR (black) and Raman (red) spectra of CoPc calculated using (a) B3LYP, shown in comparison to (b) experimental spectra represented using arbitrary intensity units. The IR spectrum is reproduced from ref 80, while the SERS data is taken from ref 81. (c) The B3LYP frequencies are also shown as bar graphs and broadened spectrum in comparison to (d) DFT+U results. The normal modes calculated with Ueff = 4 eV are shown in yellow, while the Ueff = 2 eV results are shaded in gray. Only the frequency window 500−1800 cm−1 is shown.

Ueff = 3 eV (0.59 μB) match the values determined by B3LYP (0.57 μB) and, respectively, HSE06 (0.59 μB). For FePc and CoPc, the magnetic moments per metal atom, as calculated by B3LYP and HSE06, are close in value to each other and larger than any of the GGA+U results. In FePc, the HSE06 calculated value of 2.12 μB (2.10 μB by B3LYP) is relatively higher than the moment calculated using Ueff = 5 eV (2.03 μB). Similarly, for CoPc, the 1.10 μB determined by HSE06 (1.09 μB by B3LYP) is slightly larger than the maximum of 1.06 μB determined by Ueff = 5 eV. Interestingly, a general trend in all four TMPcs is that HSE06 predicts slightly larger magnetic moments, compared to B3LYP. However, the difference is visibly larger in CuPc and, especially, in MnPc, reinforcing the conclusion that the manner in which exchange is included plays a more important role on the magnetic moments per metal atom for CuPc and MnPc than for CoPc and FePc. 3.2. Structural Properties. Table 1 contains the calculated M−Niso bond lengths obtained by GGA+U using different values of Ueff, in comparison to hybrid functional results and experimental values. First of all, what becomes evident from the table is that the D4h symmetry of MnPc is reduced to D2h by the hybrid calculations and the x- and y-axes are no longer equivalent. In the case of FePc, the M−Niso bonds along the xaxis and, respectively, y-axis differ by