DFT approaches to transport calculations in magnetic single-molecule devices

Electron transport properties of single-molecule devices based on the [Fe(tzpy)2(NCS)2] complex placed between two gold electrodes have been explored using three different atomistic DFT methods. This kind of single-molecule devices is quite appealing because they can present magnetoresistance effects at room temperature. The three employed computational approaches are: (i) self-consistent non-equilibrium Green functions (NEGF) with periodic models that can be described as the most accurate between the state-of-art methods, and two non-self-consistent NEGF approaches using either periodic or non-periodic description of the electrodes (ii and iii). The analysis of the transmission spectra obtained with the three methods indicates that they provide similar qualitative results. To obtain a reasonable agreement with the experimental data, it is mandatory to employ density functionals beyond the commonly employed GGA (i.e., hybrid functionals) or to include on-site corrections for the Coulomb repulsion (GGA+U method).


Introduction
The field of Molecular Electronics has been developed with the goal to provide a miniaturization of the electronic devices beyond the limit of the silicon technology. [1] Nowadays, the development of new devices that can overcome the Moore law is mandatory because the actual size of the gates (14 nm in last technologies) of the silicon-based transistors is close to the limits. Thus, new technologies "More than Moore" and "beyond CMOS" [2] are required to achieve future generations of computing devices. Under this perspective, heterogeneous devices replacing some silicon components by nanometric 2D or 1D-systems and molecules are a good alternative.
Furthermore, the possibilities for controlling electron spin in molecular systems open many opportunities and challenges in the Molecular Spintronics field. [3,4] Several examples of devices based on multilayer systems exploiting spin properties have been developed (e.g. spin filter, spin valves, negative differential resistance devices). [5][6][7] Most of the devices are based on a nonmagnetic layer sandwiched between two magnetic electrodes showing among other magnetoresistance properties. Such property consists in a change of the conductivity when the magnetic polarization of one of the magnetic electrodes in switched. The use of molecules in this research field is relatively new and most of the experiments were performed far from room conditions, such as ultrahigh vacuum and low temperatures. [8][9][10][11] The development of new devices has been done basically using break-junction or scanning tunneling microscope (STM) experiments with the magnetic molecule (mainly single-molecule magnets and spin-crossover systems) placed between the two metal electrodes. Recently, some of us reported molecular-based devices displaying room-temperature spin-dependent transport; single-molecule STM junctions built by bridging individual small magnetic molecules such as spin-crossover Fe II complexes [Fe(tzpy) 2 (NCX) 2 ], (X = S or Se) in a high-spin S = 2 configuration deposited on a gold substrate using a magnetic nickel tip. [12] From the computational point of view, many theoretical models have been proposed for the study of Molecular Electronic systems. Atomistic calculations remain restricted to DFT methods, [13][14][15] due to the complexity and large number of atoms (two metal electrodes plus the magnetic molecule) present in the devices. [16,17] DFT approaches are commonly based on the assumption that the electrons travel through the molecule rapidly (tunneling) without causing reduction or oxidation of the molecule. In the experiments, such conditions can be achieved with small molecules and good molecule-electrode contacts. The opposite mechanism is the Coulomb blockade regime where the conduction electron remains captured in the molecule (usually with larger molecules and bad contacts) has been mainly explored using Master Equation formalism but also time-dependent DFT methods. [18,14] Furthermore, DFT-based approaches consider that the electron transport in the scattering region (molecule and surface electrode in the contacts) is coherent (no change in the wavefunction phase) and elastic if not change in the energy in the transport process despite that some inelastic corrections can be mainly included to take into account the electron-phonon interactions.
The most common procedure to analyze the phase coherent transport in molecular transport junctions is based on the work of Landauer, Imry and Buttiker. [19] Their expression for the conductance (G) for a given system with current I and voltage V is, where T ii , e and h are the transmission through the channel i, the electron charge and the Planck's constant. Usually, the pre-factor combination of constants 2e 2 /h is defined as the quantum of conductance, G 0 (conductance through a single-atom metal wire). From the practical point of view, the sum in Eq. 1 is considered in terms of the molecular orbitals of the molecule that can provide an electron pathway channel between the two electrodes. The energy of such orbitals must be closer to the Fermi level of the metal electrodes than the applied voltage. During the travel through the molecule, there is a scattering process but there are not significant changes in the electronic structure of the molecule (tunneling) and the process is energy conserving (in the original Landauer formulation, the scattering was always considered elastic). Eq. 1 cannot be directly applied by quantum chemistry methods. However, using non-equilibrium Green functions (NEGF), [20] the current (I) can be expressed as the integral over some voltage-and energy-dependent magnitudes.
( * /0 +1 Γ 3 0, 5 6 7 0, 5 Γ 8 0, 5 6 9 0, 5 (; 3 0, 5 − ; 8 0, 5 ) (2) where G r (G a ) is the retarded (advanced) Green function, Γ is the spectral density of the electrodes (twice of the imaginary component of the self-energies) and f is the Fermi distribution of the electrodes (left and right). From Eq. 2, the transmission can be expressed as: T E, V = +1 Γ 3 0, 5 6 7 0, 5 Γ 8 0, 5 6 9 0, 5 The terms needed to apply Eq. 3 can be extracted from the corresponding Hamiltonian, and the overlap matrices from any electronic structure method, for instance, extended-Hückel (tightbinding) or DFT. Thus, the retarded Green's function matrix for the whole system:    The goal of this paper is to analyze three different DFT approaches and their accuracy to study coherent elastic spin-polarized transport through a magnetic single-molecule device (see Figure 1, the spin-crossover Fe II complex, [Fe(tzpy) 2 (NCS) 2 ], mentioned above between gold electrodes): (i) The first approach employs a combination of NEGF technique with density functional theory to obtain the self-consistent mean-field Hamiltonian of the system subject to a finite bias voltage and from it, the Green's functions providing the non-equilibrium electronic density and current. [21,14] This approach is available in several computer codes, for instance Transiesta, [21] ADF-BAND, [22,23] Smeagol [24] or ATK, [25] and allows a periodic description of the electrodes, normally associated with large computational demands. Magnetic molecules can pose additional complications for this approach, as correct orbital occupations must be retained in calculations see ref. [27]) will result in small energy gaps between empty and occupied orbitals, thus, giving wrong transmission curves and conductance values that are often too large. This can be improved by using more sophisticate functionals that partially remove such error (for instance, hybrid or long-range corrected functionals). (ii) The second approach tries to solve the main drawbacks of the first approach but paying the price of describing the electrodes as a relatively small finite cluster of metal atoms, which allows to perform a common quantum chemical molecular calculation on them (Gaussian, [28] Q-Chem, [29] ADF, [30] Turbomole [31] or FHIaims [32]). Thus, hybrid or long-range corrected exchange-correlation functionals can be employed to provide better orbital energies, and consequently transport properties. This kind of computer codes have also control over the electronic structure of the magnetic molecule once placed in between the electrodes to have correct open shell electron distribution. The Green's functions are calculated in the less-accurate wide band limit approximation, which assumes that the density of states (DOS) of the electrode is independent of the energy. [22] Thus, the Green's functions of the electrode surface (in Eq. 6) are just a constant value under the Fermi level. In these codes, the transport module is often implemented as a post-processing tool (Artaios, [33] AITRANSS [34]) and are usually restricted to zero-bias, however, other codes as Alacant [35] also allow a fully self-consistent NEGF approach implemented in the Gaussian code. (iii) Finally, the third approach adopted in the GOLLUM computer code [36] that uses periodic tight-binding calculations or mapping DFT Hamiltonian matrix (calculated with other code, for instance Siesta [37,38] that can incorporate DFT+U and dispersion corrected functionals and also can provide a good control of the electronic structure of magnetic open-shell systems) with a tight-binding approach. Next in a postprocessing step with simpler equilibrium transport theory can provide a large variety of transport properties with less computational resources than with the self-consistent NEGF codes.

Computational details
Calculations were performed with ATK, [25] Gaussian [28]-Artaios [33] and Siesta [38]-Gollum [36] packages. The self-consistent NEGF were performed with the ATK code using periodic models using PBE [39] and PBE+U functionals and double-zeta basis set with polarization for all the elements with the exception of gold that a single-zeta basis with polarization was used combined with an 11-electron pseudopotential. The periodic model structure has the thiocyanato S atoms located in threefold hollow sites of the Au(111) surface [40] (distance S-surface of 2.5 Å, Figure 1).
The calculations were carried out with the molecule sandwiched between three gold layers with a 6 × 6 surface unit cell.
The Gaussian calculations were performed with the release D01 using the pure GGA (PBE), [39] hybrid (B3LYP), [41] hybrid meta-GGA (TPSSh), [42] and long-range corrected (ωB97X) [43] exchange-correlation functionals with a lanl2dz basis set including pseudopotentials for all the elements. [44] The choice of the TPSSh functional is due to the well-know problems of most of DFT functionals [45] to reproduce correctly the energy difference of high-and low-spin states of spin-crossover complexes and such exchange-correlation functional provides the best results. [46][47][48] The post processing procedure to obtain the transmission curves was performed with the Artaios code (version 1.9). [33] The non periodic model has the same structure that the periodic model used with ATK but the gold electrode is cut just keeping 22 gold atoms (3,7,12 atoms in the first, second and third layer, respectively).
Calculations with the Siesta code [38,37] were performed using the version 3.0 using PBE [39] and PBE+U functionals with and double-zeta basis set for all the elements with the exception of gold that a single-zeta basis was used combined with an 1-electron pseudopotential.

Results and Discussion
Previously, some of us have analyzed using PBE functional and periodic calculations with the selfconsistent NEGF Transiesta code the transport properties of the high and low-spin states of the [Fe(tzpy) 2 (NCS) 2 ] complex placed between two gold electrodes. [40] The high-spin state (S = 2, t 2g 4 e g 2 ), with five alpha and one beta electrons, shows a much higher conductance than the low spin (S = 0, t 2g 6 e g 0 ). [40,50] The reason for the high-spin large conductivity, considering that there is only one beta electron in a t 2g orbital, can be found in the fact that the orbital bearing this electron, as well as the first unoccupied beta t 2g orbital, are those closer in energy to the Fermi level of the gold electrodes providing effective channels for transport through the molecule (see Figure 2).
However, in the low-spin state, both sets of orbitals (t 2g and e g ) are relatively distant from the Fermi level, thus resulting in a lower conductivity. As in the high-spin state, the transport is due to the beta orbitals. The current is highly spin polarized of minority beta carriers being such property crucial for the magnetoresistance behavior found experimentally using STM measurements of such systems. [12] Self-consistent NEGF with periodic models The calculations were performed using the ATK code [25] and the PBE [39] (and PBE+U with U=4.0 eV a recommended value for Fe II complexes [51]) functionals. The transmission spectra for the high-spin S=2 Fe II complex for PBE and PBE+U functionals are represented in Fig. 2. There is a clear influence of the inclusion of the electron repulsion through the U parameter resulting in a much larger gap between occupied and empty orbitals. Thus, the t 2g levels are more distant of the Fermi level of the electrodes and also logically a smaller intensity is obtained for such method (see     (blue) is for rightward transport; π (red) is for leftward transport; and π/2 is transport perpendicular to the electrode plane.

Non Self-Consistent NEGF with cluster models
As mentioned above, the non self-consistent NEGF methods using cluster structural methods have the limitations of the structural model and also the simplifications in the calculation of the surface Green functions. However, it is possible to employ common Quantum Chemistry codes (i.e., Gaussian09 [28]) opening the possibility to use a wide range of exchange-correlation functionals and to extract the transport properties (only transmission spectra) in a post-processing procedure (Artaios code [33]). In the previous section, we have shown the usual overestimation of the transport properties provided by the GGA functionals. Thus, in our analysis we will compare four different functionals: (i) PBE functional Perdew, 1996 #51} to have the same reference than in the self-consistent NEGF periodic method, (ii) the popular hybrid B3LYP functional, [41] (iii) the hybrid meta-GGA functional, TPSSh, [42] usually recommended to provide a good description of spin state stabilities in metal complexes and (vi) the long-range corrected ωB97X functional [43] that should improve the asymptotic behavior and consequently, to reduce the self-interaction error to give a better description of the excitations.
The comparison of the PBE transmission spectra (Figures 2 and 6) shows very similar results with two broad beta transmission peaks close to the Fermi level (indicating strong contact interaction with the gold electrode) with a distance between them around 0.5 eV (2 eV for PBE+U U=4.0 eV) while the alpha gap is around 2 eV. The position of the beta transmission peaks is the key parameter for the transport properties; thus, the gap between both peaks (with large contribution of t 2g -like metal orbitals, see Figure 6) is around 3, 2.3 and 8 eV, for B3LYP, TPSSh and ωB97X functionals. Thus, the two hybrid functionals give similar results (slightly higher energy gap for B3LYP due to the larger exact-type exchange contribution) and close to those obtained with the self-consistent NEGF periodic PBE+U calculations. However, the long-range corrected ωB97X functional gives a much higher energy gap. Concerning the position of these two peaks compared with the Fermi level, the empty t 2g -like orbital remains relatively close with all the functionals while the occupied one is considerably shifted to low energy values. Thus, such result would indicate that the transport would be through the unoccupied levels as it is usually found for this kind of single-molecule devices. properties are also plotted.

Non Self-Consistent NEGF with periodic models
This method allows to perform the calculation of transport properties in a post-processing way but using periodic boundary conditions for the description of the electrodes. Thus, the option to improve the results with DFT+U methods is available using regular calculations with the Siesta code [38,37] and obtaining lately the transport properties with the Gollum package. [36] Furthermore, the fact to be based on a non self-consistent approach results in a considerable  It is worth noting that another representation of transport properties commonly employed with the Gollum program is to analyze the dependence of the G conductance with the Fermi level value. [52,53] This representation (see Fig. 8) is often employed because it is well-known drawback that the calculated Fermi level values with common DFT functionals are usually non-accurate, especially taking into account that the transport properties are highly dependent of a small energy shift of the transmission peaks. Thus, it is possible to easily choose the Fermi level value that provides a reasonable agreement with the experimental data if not the Fermi level is fixed at zero.
The I/V characteristics can also be calculated (see Fig. 9) with such post-processing approach leading to almost identical results than those obtained with the self-consistent NEGF approach.

Concluding Remarks
Theoretical methods can provide very useful information about the transport mechanism in singlemolecule devices. The presence of magnetic molecules in such devices opens the field of Molecular Spintronics to new physical properties. We have analyzed three different theoretical approaches to study devices based on the [Fe(tzpy) 2 (NCS) 2 ] complex that shows room temperature magnetoresistance in single-molecule devices. Due to the small size of the molecule and the relative good molecule-gold electrode contact, the transport mechanism in such system is basically by tunneling. Hence, such mechanism can be analyzed using theoretical methods using atomistic DFT calculations using the non-equilibrium Green function approach to calculate transmission and current properties using Landauer equation. The most accurate approach is the self-consistent NEGF using a periodic description of the electrodes. The calculations were performed using PBE and PBE+U (U=4.0 eV for Fe d orbitals). Due to the well-know drawback of the large selfinteraction error of the PBE functional, the presence of occupied and empty molecular levels too close to the Fermi level of the electrode results in a considerable overestimation of the current.
Thus, the PBE+U approach improves considerably the results in such methodology. The inclusion of the +U correction shifts the orbitals far of the Fermi level, increasing the energy gap, and consequently the overestimation of the current obtained with the PBE functional is corrected.
The second employed theoretical approach allows to perform the calculations with a general computer code where many different flavors of the functional are available (i.e. Gaussian09). The transport properties are non self-consistently calculated using NEGF with a post-processing tool using small finite metal clusters to model the electrodes. The results using the same functional that the first approach provide similar results and again, the use of hybrid functionals to reduce the selfinteraction error is needed to obtain a similar transmission curves that the obtained with the periodic self-consistent NEGF DFT+U calculations.
Finally, the third approach is also based on a periodic description of the electrodes but also with a post-processing tool the transport properties are obtained in a non self-consistent NEGF method.
The comparison of the same functionals (PBE and PBE+U) than with the self-consistent NEGF periodic method give actually very similar transmission curves with considerable smaller computational resources. Hence, we can conclude that qualitatively the three theoretical approaches can achieve a right semi-quantitative description if the appropriate functional is employed and, qualitatively even the simplest GGA functional provides a correct description of the orbitals involved in the transport channels despite a strong overestimation of the calculated current.