Platinum atomic contacts: from tunneling to contact

TL;DR

This paper provides a detailed theoretical analysis of electronic transport in platinum nanocontacts, emphasizing the importance of technical details like system size and basis set choice in accurately modeling tunneling and contact regimes.

Contribution

It highlights the critical influence of system size and basis set selection on tunneling conductance calculations in platinum nanocontacts, offering practical recommendations for accurate modeling.

Findings

Proper system size prevents unphysical charge oscillations.

Large-cutoff-radius basis sets improve conductance decay modeling.

Technical details have minimal impact in the contact regime.

Abstract

We present a theoretical study of the electronic transport through Pt nanocontacts. We show that the analysis of the tunnelling regime requires a very careful treatment of the technical details. For instance, an insufficient size of the system can cause unphysical charge oscillations to arise along the transport direction; moreover, the use of an inappropriate basis set can deviate the distance dependence of the conductance from the expected exponential trend. While the conductance decay can be either corrected by employing ghost atoms or a large-cutoff-radius basis set, the same does not apply to the corrugation, for which only the second option is recommended. Interestingly, these details were not found to have a remarkable impact in the contact regime. These findings are important for theoretical studies of distance-dependent phenomena in scanning-probe and break-junction experiments.

Figures

e and Fig.8 of ref Kizuka and Monna ([2009](#bib.bib61))), which is closer to the LANl2dz value.

Equations

$$T_{\sigma }(E)=\frac{1}{V_c}{\int }_{BZ}d{k}^3{T}_{\sigma }^k(E)$$

$$T_{\sigma }^{(k)}(E)=Tr[{\Gamma }_{\sigma ,L}^{(k)}(E)G_{\sigma ,C}^{(k)}(E+iϵ){\Gamma }_{\sigma ,R}^{(k)}(E)G_{\sigma ,C}^{(k)}(E-iϵ)]$$

$$G_{\sigma ,C}^{(k)}(Z)=[Z{S}_C^{(k)}-H_{\sigma ,C}^{(k)}-{\Sigma }_{\sigma ,L}^{(k)}(Z)-{\Sigma }_{\sigma ,R}^{(k)}(Z){]}^{-1}$$

$${\Sigma }_{\sigma ,L}^k(Z)=(Z{S}_{CL1}^{(k)}-H_{\sigma ,CL1}^{(k)})G_{\sigma ,L(Z)}^{(k)}(Z{S}_{L1C}^{(k)}-H_{\sigma ,L1C}^{(k)})$$

Full text

I Introduction

Over the past decade, metal nanocontacts have attracted remarkable interest as they make it possible to study the electronic, mechanical and thermal properties of metal atoms in low-coordination conditions, Sattler (2010); Evangeli et al. (2015); Cuevas and Scheer (2010); Jelínek et al. (2005, 2008); Evangeli et al. (2015); Müller et al. (2016) besides offering the possibility of being operated as two-terminal switches.Schirm et al. (2013) Furthermore, a thorough understanding of their electron transport properties is crucial in molecular electronics, where the total transmission through metal-molecule-metal junctions often includes contribution from metal-to-metal current. Cuevas and Scheer (2010) Investigating metal contacts at the nanoscale has been possible especially thanks to the refinement of scanning tunnelling microscopy (STM) and mechanically-controlled-break-junction (MCBJ) techniques. In particular, studies on Pt contacts have revealed interesting effects. For example, it was observed that, upon elongation, the conductance oscillates and the number of scattering channels changes Fernández-Rossier et al. (2005); García-Suárez et al. (2005); Vardimon et al. (2014). In addition, it was suggested that Pt is a good candidate as electrode material in molecular junctions as it makes rectification independent of stretching distances Parashar et al. (2014). The trapping of molecules in Pt-Pt junctions, such as CO, has also been studied quite extensively and is still the object of interesting analyses Makk et al. (2012). Nevertheless, the amount of theoretical studies on electron transport through Pt nanocontacts García et al. (2004); Strange et al. (2006, 2008); Cuevas et al. (2003); García-Suárez et al. (2005); Zhang et al. (2010); Nielsen et al. (2003, 2002); Wu et al. (2014); Thygesen and Jacobsen (2005); Pauly et al. (2006); Smogunov et al. (2008); Ono (2009) still remains limited compared, for instance, to gold. The latter is more widely used in experiments due to its noble character and, consequently, has often been employed as staple electrode material in theoretical simulations. Furthermore, the Au electronic structure is particularly convenient for transport calculations, as we proceed to explain. The most common tool to compute electron transmission is the combined DFT(density functional theory)+NEGF(Non-Equilibrium Green’s Functions) method, Pauly F and G (2008); Palacios et al. ; Brandbyge et al. (2002); Rocha et al. (2006); Ozaki et al. (2010); Xue et al. (2001); Ferrer et al. (2014); Garcia-Lekue et al. (2015) which is by now implemented in several codes. Many of them include self consistent cycles in which convergence must be reached for input and output charge density. Such a task is easier on Au structures than on their Pt analogues, due to the fact that the Au valence band consists of the 6$s$ orbital mainly, while that of Pt consists of both the 6$s$ and 5$d$ bands. The limitations and inaccuracies of DFT are well known, as is the huge effort which is currently being made to try to overcome them by applying corrections or using alternative theoretical techniques.Quek et al. (2007) Nevertheless, the DFT+NEGF method is still widely used, especially for large systems. In this piece of work, we will show theoretical results based on benchmark electron transport calculations on Pt nanocontacts. We will discuss how to perform them in an accurate manner; furthermore, we will examine whether strategic approximations can be made to make these calculations less cumbersome without jeopardizing the quality of the final result. In particular, we will focus on the transition from tunnelling to contact regime, which has not been studied to the same extent as contacted systems Limot et al. (2005); Garcia-Lekue and Wang (2010) but which is extremely important when studying distance-dependent phenomena in STM or MCBJ experiments. For instance, it was shown that a change in conductance due to inelastic effects makes it possible to characterize the crossover from tunneling to contact; Frederiksen et al. (2007) it was also observed, in spin-polarized system, that the transmission probabilities of the eigenchannels present a nonmonotonic behavior as a function of the tip-adatom separation. Polok et al. (2011)

II Methodology

We performed electron transport calculations by means of the NEGF extension of the OpenMX code Ozaki et al. (2010), which is by now widely used. Sharma et al. (2013); Okuno and Ozaki (2012); Ansari et al. (2012); Lan et al. (2015); Hashmi et al. (2015); Jippo et al. (2015) The model implemented in this software makes use of periodic boundary conditions. The basis sets employed consist of linear combinations of atomic orbitals generated by a confinement scheme which yields wave functions with zero value beyond a chosen cutoff radius. Ozaki (2003); Ozaki and Kino (2004) For all elements, basis sets which have already been conveniently predefined and optimized are provided by the software distribution. Following the OpenMX notation, we will define the basis sets by combinations of $l{x}$ sets, $l$ being the orbital and $x$ being the number of primitive orbitals used for the construction.

To describe the Pt electronic structure we either used the minimal basis set s1p1d1 (containing 5${p}$, 5${d}$ and 6${s}$ orbitals) or the s2p2d1 basis set with a cutoff radius of 7 Bohr, if not stated otherwise. Notice that, by constructing the $p$ orbital with two primitive functions, we have taken polarization into account, since 6${p}$ orbital is the lowest unoccupied orbital. In all cases, the exchange-correlation functional GGA (PBE13) Perdew et al. (1997) was employed and the integration over $k$ space was numerically performed by using a Monkhorst-Pack mesh Monkhorst and Pack (1976), where the number of $k$ points was varied depending on the specific case. For the numerical integrations and the solution of the Poisson equation an energy cutoff of 200 Ry was used. All calculations were converged until the deviation of the energy eigenvalue reached a value lower than 10*-6* H.

III Results and Discussion

In OpenMX, as in most DFT+NEGF based codes, the whole structure is thought of as divided into three regions, namely left (L), central (C) and right (R) (Fig.1). The electronic transmission is evaluated by the Landauer formula. Defining k as the Bloch wave vector and $\sigma$ as the spin index, the spin resolved transmission is given by

[TABLE]

where $T_{\sigma}^{k}$ is the $k$ resolved transmission defined by

[TABLE]

$\Gamma_{\sigma,L}$ and $\Gamma_{\sigma,R}$ are the coupling matrices to the left and right lead. The Green’s function of the central C region can be written by

[TABLE]

where the self-energies $\Sigma^{(k)}_{\sigma,L}$ and $\Sigma^{(k)}_{\sigma,R}$ are defined by

[TABLE]

More details about these equations can be found in Ref. Ozaki et al., 2010. $G^{(k)}_{\sigma,L(Z)}$ is the surface Green function, which is calculated by a preliminary standard DFT calculation. This is performed on a bulk structure which has the same chemical composition of the electrode and which is repeated periodically in all three space directions.

The first technical issue that has to be addressed is how large each region has to be along the transport direction in order to get a converged value of the conductance. This question is common to other NEGF codes such as tranSIESTA Brandbyge et al. (2002) and ATK ATK , where convergence has been studied as a function of the size of the central area. Due to differences in the generation of the basis sets between these codes and OpenMX, such as in the definition of the core potential in the confinement scheme Junquera et al. (2001); Ozaki (2003); Ozaki and Kino (2004), it is by no means obvious that the very same criteria should apply to our systems. Therefore, we have applied systematic changes to the size of both the leads and the central part of a Pt$-$Pt junction and checked how these changes affect to the electronic structure, focusing in particular to the presence of unphysical charge oscillations that can arise close to the leads and which can ultimately affect the computed conductance.

As a first example, we shall analyze the case of the transport through two Pt slabs separated by 5 Å  distance along the $\left\{111\right\}$ direction. A periodically repeated 1$\times$1 unit cell was considered, with an 8$\times$8$\times$1 $k$-point grid (Fig.1). Notice that the coordination of all atoms is large enough for magnetism to be absent Fernández-Rossier et al. (2005). Fig. 2 shows the Mulliken population on each atom along the transport direction for different combinations $l-c-r$ for the number of layers in the central region ($C$) and in the left ($L$) and right leads ($R$) (red numbers in panels a-e). More precisely, $l$ and $r$ indicate the number of layers in the bulk structure, employed for the preliminary left and right electrode calculation, which are repeated infinitely along the transport direction. In Fig. 2, the three areas are separated by vertical red lines, while the dashed green line indicates the vacuum separation. The Mulliken charges are shown for both s1p1d1 (left) and s2p2d1 (right) basis set. The avarage value of 16 (higher than the 10-electron valence of Pt) is due to having included the electrons of the 5$p$ orbital. It can be observed that strong charge oscillations (up to 0.8 $e$) appear in proximity of the leads, especially with the s2p2d1 basis set.This suggests that, while an improvement of the basis set is naively expected to yield better-quality results, in this case it can actually worsen other details such as the charge population (it is worth noting that improving the quality of the basis set was found to give qualitatively wrong conclusions also in previous NEGF work for other reasons) Herrmann et al. (2010). The presence of the fluctuations does not seem to be related to the separation distance between the leads, as it was also found to take place, for instance, in the contact regime. Similar oscillations have also been observed in NEGF calculations for carbon chains (Fig.3 of Ref. Ozaki et al., 2010), albeit with smaller amplitude. We also found that replacing Pt with Au causes a reduction of the oscillations (see the corresponding Mulliken population for the system $3-6-3$ in the inset of Fig. 2). In this case, the population did not seem to be affected by an improvement of the basis set as much as for Pt, probably because of the different electronic structure. Fig. 2 shows that, in order to keep a significant portion of the central region unaffected by such oscillations, it is necessary to either increase $l$ ($r$)( as in panel d-e) or $c$ (as in panel b-c). However, the second option is more recommendable as the first one causes the computational time of the NEGF run to increase considerably.

Subsequently, we proceeded to check how these charge oscillations affect the energy dependence of the electronic transmission. Fig. 3 shows the transmission curves for the systems (a)-(e) of the right column of Fig. 2, calculated with the s2p2d1 basis set and a 48$\times$48 $k$-point grid at three different separation distances (0, 2 and 5 Å). It can be observed how, at 0 and 2 Å, no remarkable differences appear in the energy range around the Fermi level. At 5 Å instead, in the same region the transmission for the 3-6-3 system is strikingly different from the others, as the transmission in the energy region below is generally lower than in the curves of the other $l-c-r$ sets. Although the low$-$bias conductance values of the five transmission curves are not very different from each other, overall the results we presented so far indicate that the 3-6-3 system is not large enough and either the size of the leads or of the central region must be extended as in the other cases. Since extending the leads was found to increase the computational time more than increasing the central part did, we propose the $3-12-3$ structure as a good compromise.

We now turn to analyze the distance dependence of the conductance and how it is affected by differences in the basis sets. In the inset of Fig. 4, we show the geometry we employed. We considered two cases, in which the apex of the tip of the upper lead falls on an hcp and on a top site of the lower lead. The structure is a 4$\times$4 periodic cell and the central area comprises 6 layers on each side, while each lead consists of 3 layers. In both cases, we increased the separation distance from the initial value of 2.6 Å and calculated the conductance stepwise. In Fig.  4, we compare the natural logarithm ($ln$) of the conductance values obtained in both cases. We adopted three different cutoff radii for the basis set confinement (7, 9 and 11 Bohr). In the intermediate range 3.6 - 5.6 Å, all curves approximately present a linear shape (corresponding to the expected exponential decay of the conductance) and no significative quantitative differences are visible. At closer distances we can observe the typical deviation from a linear behaviour, which is expected in the contact regime Garcia-Lekue and Wang (2010). For the shortest cutoff radius of 7 Bohr, the curves show deviation from linearity also beyond 5.6 Å. This problem is well known and is usually solved by inserting ghost atoms in the vacuum region (which can occasionally hamper the scf convergence) or by using plane-wave-based methods Garcia-Lekue and Wang (2010); García-Gil et al. (2009); Garcia-Lekue et al. (2015). A comparison between the use of ghost atoms and increasing cutoff radius was previously performed by Siesta calculations García-Gil et al. (2009) for what concerns the spatial behaviour of silver wave functions. There, the insertion of ghost atoms was claimed to provide closer results to those obtained with plane waves calculations. However, instead of the wave functions, we chose to focus the comparison on a physical quantities such as the conductance and the corrugation. Notice also that, despite the thorough study of Ref. García-Gil et al., 2009, the results there obtained cannot be automatically applied to our system because of the differences in the confinement scheme, as mentioned above. In Fig. 4, we also show the results obtained by inserting ghost atoms. It can be noticed that the ensuing conductance values are very similar to those obtained by increasing the cutoff radii to 11 Bohr. Moreover, the employment of ghost atoms did not show any particular advantage concerning, for instance, the computational time needed to achieve convergence. This indicates that increasing the cutoff radius can indeed be a viable alternative solution to this well-known problem, as it yields reliable results without the insertion of additional elements in the geometry. We will also show in the following that the employment of ghost atoms can lead to a wrong evaluation of the corrugation in the tunneling regime. In Fig. 5, we compare, for each cutoff radius and for the ghost-atom case, the conductance values obtained in the top geometry with that obtained in the hollow at the same separation distance. Interpolation curves obtained by the Akima spline method in the range 2.6-4.6 Å are also shown (dashed lines). In the inset of the bottom right panel of the same figure we report the corrugation as a function of distance for all four sets. This quantity was calculated from the interpolation curves by evaluating, at each conductance value, the difference in tip-surface distance between the corresponding top and hollow geometries. Such a difference (the corrugation) has been plotted against the average tip-surface distance. In all four cases, the corrugation is negative and follows a linear trend up to about a distance of 3.5 Å. Beyond this value, however, while the curves for cutoff radii of 7, 9 and 11 Bohr reach positive values and saturate at 0.1 Å approximately, the corrugation relative to the insertion of ghost atoms fails to reproduce this inversion, remaining negative. The origin of such a discrepancy can be spotted in the main four panels of Fig. 5: in the ghost case, the conductance for the hollow position is always higher than for top, whereas, in the other three cases, an inversion takes place at 3.5 Å. This change at intermediate distances is indeed expected, as it has also been observed in other systems Kim and Hasegawa (2015).

Finally, in order to assess the reliability of our results, we explored how changes in the number of primitive functions employed to construct the basis sets affect the conductance and how these results compare with those obtained by other models. To this aim, we chose the same geometry (depicted in the inset of Fig. 6) as in Ref. Thygesen and Jacobsen, 2005, where the Pt electronic structure was described by means of Wannier functions. In our system, however, we added three more layers on each side of the central part in order to avoid the charge oscillations discussed above (Fig. 2). The distance between the two tip apices is 2.8 Å, hence in the regime in which the conductance is not affected by the cutoff radius. Note that, for a similar structure, the lower coordination of the tip apex was found to give rise to spin polarization. Korytár and Evers (2013) However, for the present purpose we chose to neglect it since including it would worsen the numerical performance considerably and would not change our main conclusions. We either used the s1p1d1 or s2p2d1 basis set for all atoms. However, we also considered a case in which not all atoms were treated in the same way: there, the atoms in the tip and the innermost layer on each side were described by either the s2p2d1 or s2p2d2 set, whereas the s1p1d1 set was used for the rest. This mixed strategy is quite frequent in theoretical studies of metal nanocontacts as it allows a good description of the central relevant region at a reasonable computational cost Wu et al. (2014); Zotti et al. (2011). For all four cases, the self consistent cycle was performed with a 4$\times$4 $k$-point grid; subsequently, the transmission, reported in Fig. 6 as a function of energy, was calculated using the same 6$\times$6 $k$-point sampling as in Ref. Thygesen and Jacobsen, 2005. The four curves (see Fig. 6) show, right below the Fermi level, two peaks stemming from the expected contribution of the $s$ and $d$ states Cuevas et al. (2003) and from the loss of degeneracy of the d states due to symmetry reasons Nielsen et al. (2002). In the critical region around the Fermi level, the s1p1d1 curve shows a slight downshift (0.05 eV) with respect to the other cases, in agreement with what was observed in Fig. 3 of Ref. Strange et al., 2008. In this range all curves appear similar to those reported in Ref. Thygesen and Jacobsen, 2005, although shifted down in energy by around 0.25 eV. In Appendix A, further comparisons with other models are shown. Overall, we can conclude that, at the Fermi level and in the contact regime, detail differences in the basis sets do not affect either the general shape of the transmission curve or the energetic orbital alignment considerably. More noticeable changes were observed, instead, at higher energies, which however do not affect the low bias conductance.

IV Conclusions

We studied the electronic transport through Pt nanocontacts in the $\left\{111\right\}$ orientation by using the combined DFT+NEGF methods as implemented in OpenMX. We showed that, while the analysis of the contact regime is not remarkably affected by details in the leads and in the construction of the basis sets, this does not hold for the tunneling regime in many ways. For instance, at large separation distances between the two electrodes forming the contact, it is important to extend either the central or the lead region to a number of layers so as to avoid the formation of unphysical charge oscillations in the region of interest. We also showed that, in this regime, it is necessary to use basis sets with a large cutoff radius to describe both the conductance decay and the corrugation accurately. While the former could also be reproduced by employing ghost atoms, the same strategy does not apply to the latter, for which ghost atoms fail to reproduce the inversion in the corrugation expected at intermediate distances between top and hollow position. These findings suggest that similar preliminary checks should be made for other materials, whenever the performance of NEGF calculations aims at the study of distance-dependent effects. This concerns, for example, the analysis of scanning-probe experiments. Ultimately, we proved the robustness of our results by comparing them with those obtained by other NEGF-based models.

V Acknowledgments

We thank Eduardo Hernandez, Juan Carlos Cuevas and Juan Jose Palacios for fruitful discussions and M. Todorovic for help with the DFT calculations. LAZ was funded by the Spanish MINECO through the grant MAT2014-58982-JIN. RP acknowledges financial support from MINECO projects CSD2010-00024, MAT2011-023627, MAT2014-54484-P, and MDM-2014-0377. Computer time provided by the Spanish Supercomputing Network (RES) at the MareNostrum III (BSC) and Tirant is gratefully acknowledged.

Appendix A Comparison with ANT

To further corroborate the robustness of our results, we investigated how the energy dependence of the transmission is affected by detail differences in the DFT+NEGF model employed. To this aim, we repeated our calculations by the quantum-chemistry code ANT  Jacob and Palacios (2011); Palacios et al. which is built as an interface to Gaussian Frish et al. (2003). There, metals are described by means of isolated clusters (thus no periodic boundary conditions are applied as in OpenMX), while the basis sets consist of linear combination of gaussian functions. Combined DFT+NEGF methods are by now commonplace and have been previously compared with other levels of theory Brandbyge et al. (2002). However, to the best of our knowledge, comparisons between two codes which employ such fundamentally-different descriptions of the leads are few and far between Pauly et al. (2008); Bilan et al. (2012); Strange et al. (2008). We chose to recompute the transmission for the structure of Fig. 6. To this aim, we employed a geometry consisting of two Au20 pyramidal clusters. As for the basis set, we used a CRENBS Ross et al. (1990) basis set (which includes the same orbitals as in the s1p1d1 set of OpenMX) and a LANL2DZ set Wadt and Hay (1985).

In Fig. 7, we show a comparison between the transmission curves obtained by ANT and OpenMX. It can be observed that, in the range around the Fermi level, the shape of the OpenMX s1p1d1 curve and the ANT CRENBS curve are quite similar, as well as the energy alignment. Differences arising below this range are probably due to how the core part is treated. The LANl2dz basis set was previously claimed to give reliable results due to the fact that it does not add ghost transmissions as supposidely-higher-quality basis sets do. Herrmann et al. (2010) Interestingly, the general shape of the curve obtained by OpenMX with its minimal basis set appears to be quite similar to that calculated by ANT with the LANl2dz set. This indicates that, actually, the quality of the OpenMX minimal basis set is comparable to that of well known higher-quality basis sets. However, the LANl2dz conductance values are higher by approximately 0.5 G/G0. The computed values can hardly be directly compared with the experiments, since typical MCBJ measurements provide histograms which cannot be directly related to exact detailed geometries. Experimental values around 0.5 and 1.5-2 G/G0 have been reported  Evangeli et al. (2015); Vardimon et al. (2014); Nielsen et al. (2003) which do not necessary correspond to the geometry analyzed here. Indeed, because of the presence of the anisotropic $d$ states at the Fermi level, the conductance is expected to be strongly influence by geometrical detail differences throughout the process of formation and stretching of the junction. High-resolution transmission electron microscopy indicate, for the geometry of Fig. 6, a conductance value of 4 G/G0 (Fig.7e and Fig.8 of ref Kizuka and Monna (2009)), which is closer to the LANl2dz value.

Appendix B Scf convergence in OpenMX

In the NEGF calculations as those discussed in this work, it is necessary to achieve convergence between the input and output charge densities (let NormRD be the residual norm between these two quantities). It is well known that this is facilitated by performing a preliminary standard DFT run on the central region to obtain the corresponding charge density as a starting point. In OpenMX, convergence can be achieved by four different mixing schemes. Surprisingly, in all systems analyzed for this work, the algorithm RMM-DIIS with Kerker’s metric Kresse and Furthmuller (1996), which is slow but very efficient for standard DFT calculations, often failed to reduce NormRD in the NEGF runs. Instead, convergence was eventually reached by using the Pulay algorithm Bowler and Gillan (2000). In particular, we found it convenient to adopt a three-step procedure: we first performed the preliminary DFT calculation by the RMM-DIIS algorithm until NormRD reached values smaller than at least 0.01; we then switched to using the Pulay algorithm for a few runs (about 10); we finally switched to the NEGF technique by keeping the Pulay mixing scheme on. Interestingly, we also observed that the number of steps necessary for convergence did not always decrease by reducing NormRD to a much lower value than 0.01 (Fig.  8).