Interfacial charge density and its connection to adhesion and frictional forces
IInterfacial charge density and its connection to adhesion and frictional forces
M. Wolloch, G. Levita, P. Restuccia, and M. C. Righi
1, 2, ∗ Department of Physics, Informatics and Mathematics,University of Modena and Reggio Emilia, Via Campi, 213/A 41125 Modena, Italy CNR-Institute of Nanoscience, S3 Center, Via Campi 213/A, 41125 Modena, Italy
We derive a connection between the intrinsic tribological properties and the electronic propertiesof a solid interface. In particular, we show that the adhesion and frictional forces are dictated bythe electronic charge redistribution occurring due to the relative displacements of the two surfacesin contact. We define a figure of merit to quantify such charge redistribution and show that simplefunctional relations hold for a wide series of interactions including metallic, covalent and physicalbonds. This suggests unconventional ways of measuring friction by recording the evolution of theinterfacial electronic charge during sliding. Finally, we explain that the key mechanism to reduceadhesive friction is to inhibit the charge flow at the interface and provide examples of this mechanismin common lubricant additives.
Friction and adhesion are common phenomena that im-pact many fields from nanotechnologies to earthquakes,but their fundamental origin is still largely unknown [1].The reason resides in the fact that even for macroscopicobjects, friction and adhesion are governed by micro-scopic contacts, whererather uniformly (black solid curveof Fig. 5(d)), while for the on-top configuration signif-icantly less charge (red dashed curve of Fig. 5(d)) isconcentrated at the interface and it is less homogeneousFig. 6(b) and (d). the atomistic interactions of quantummechanical origin ultimately determine the tribologicalresponse [2]. Thus, it is of great practical and theoreti-cal importance to understand the connection between theelectronic structure and the mechanical tribological prop-erties of interfaces. At the atomic level, adhesion is dic-tated by the chemical interaction between the surfaces incontact and adhesive friction arises because this interac-tion changes as a function of the relative lateral positionof the two surfaces. In turn, the adhesion and frictionalforces can be understood by analyzing the charge density ρ in the region of the interface, and more specifically thecharge redistribution occurring when the two surfaces aremoved relative to each other, either to initially form theinterface or during sliding.Understanding the connection between the interfacialcharge density and adhesive friction is of paramount im-portance to design lubricant additives [3, 4]. However,nowadays most of the research on lubricants is conductedempirically due to the lack of predictive understand-ing, which we believe can be achieved by the analysisof the electronic interfacial properties. The functionalityof some solid and boundary lubricants is, in fact, basedon their capability of decreasing the adhesive interactionsbetween the surfaces in contact. In this Letter we showthat this functionality relies precisely on their ability toreduce the charge density at the interface.The advent of scanning probe techniques in tribology,such as friction force microscopy, has allowed scientiststo obtain friction maps between nanometer-size contactswith nano- and piconewton resolution [5]. Here we show that the friction maps directly reflect the charge densitymaps recorded during sliding. Therefore simultaneousmeasurement of the tribological and electronic interfacialproperties should be attempted.One of the most important figures of merit in tribologyis the work of separation, which corresponds to the en-ergy required to separate two surfaces from contact andis the opposite of the adhesion energy W sep = − E adh .The variation of E adh as a function of the lateral dis-placement during sliding is what causes the appearanceof frictional forces. This variation is described by a po-tential energy surface (PES), V ( x, y, z eq ), which in thedislocation community is known as the γ -surface, andwhere z eq is the surface separation at zero load. The ab-solute minimum of the PES corresponds to the adhesionenergy, E adh = V ( x eq , y eq , z eq ), while the energy differ-ence between the minimum and the maximum of the PESis referred to as the corrugation and will be denoted by∆ V in the following. This number is especially importantsince it is equivalent to the maximum amount of energyper unit area that might be dissipated by frictional pro-cesses.Historically simple sinusoidal energy profiles haveplayed a significant role in describing the elementarymechanisms of friction, where stick-slip (or continuousmotion) from minimum to minimum is analyzed [6, 7].Later, ab-initio data were used to generate these en-ergy profiles [8, 9], and most recently, the whole twodimensional PES has been used to analyze friction [10–15]. Using the whole PES allows one to identify frictionanisotropy and the easiest sliding path (or minimum en-ergy path) which carries the highest statistical weight.In 2012, Reguzzoni and coworkers used the interfacialcharge density, especially charge density difference pro-files, to gain insight into the frictional characteristics ofgraphene sliding on graphene [16]. In recent years otherpublications built on this idea [13, 17–19].In this Letter we use density functional theory (DFT)to present our discovery of a deeper connection betweenadhesion, the PES, and interfacial charge density varia- a r X i v : . [ c ond - m a t . o t h e r] J u l tions. We consider a large number of solids and find alinear relation between the amount of charge that is re-distributed during the formation of an interface, and theadhesion energy. Moreover, a simple functional relationbetween the strength of adhesion E adh and the corruga-tion of the PES ∆ V is discovered. We also explain thatone key function of lubricant materials is surface passi-vation to impede the charge flow at the interface [20–23].Finally we show that the PES corrugation and in turnadhesive friction are determined by the variation of thetotal charge density at the interface during sliding. Anexperimental verification of this observation using scan-ning probe techniques is proposed.For all calculations we used the plane wave DFT pack-age Quantum ESPRESSO [24], the computational detailscan be found in the Supplementary Material (SM) [25]. FIG. 1. Calculation of ρ diff using Fe(110) surfaces as an ex-ample. The total charge density ρ is visualized in a black-brown-white color scheme on a log10 scale. ρ diff is visualizedin a color scheme from blue (depletion) to red (accumulation)of charge. Line plots are planar averages of ρ diff and ρ respec-tively on a scale of 10 − electrons per ˚A . Fig. 1 shows how the charge displacement, or chargedensity difference ρ diff , is calculated for an example sys-tem of Fe(110) surfaces. The procedure is exactly thesame for the adhesion energy E adh , as is indicated by theformulas in Fig. 1: The charge density (energy) of thetwo parts is subtracted from the charge density (energy)of the combined interface. Total charge is plotted usinga black-brown-white color scale and ρ diff is shown on ablue-white-red scale to visualize depletion (blue) and ac-cumulation (red). Planar averages are also shown as aline profile in the direction normal to the interface allow-ing one to quickly see where most charge is redistributedupon interface formation. There the units are 10 − elec-trons per ˚A , since we divide by the surface area of thesimulation cell to compare rather uniformly (black solidcurve of Fig. 5(d)), while for the on-top configuration sig-nificantly less charge (red dashed curve of Fig. 5(d)) isconcentrated at the interface and it is less homogeneousFig. 6(b) and (d).different systems. It is immediatelyclear that the charge density of the separated slabs decaysexponentially in the vacuum. If the slabs are brought to- gether, the charge near the surfaces is depleted slightlyand accumulated in the interface region, which we defineas the space between the lowest atomic layer of the topslab (at z ) and the highest one on the bottom slab (at − z ).To quantify this redistribution of charge density we in-tegrate the absolute value of the profile of ρ diff in theinterface region, and normalize with respect to its width.We call this figure of merit, which measures both de-pletion and accumulation of charge density within theinterface ρ redist ,rather uniformly (black solid curve ofFig. 5(d)), while for the on-top configuration significantlyless charge (red dashed curve of Fig. 5(d)) is concentratedat the interface and it is less homogeneous Fig. 6(b) and(d). ρ redist = 12 z (cid:90) z − z | ρ diff | dz . (1)This redistribution is generally more important than the(usually very small) net flow of charge into the inter-face region. ρ redist corresponds to the shaded area in theline profile of Fig. 1, normalized to the interface width2 z . Because of this normalization, which is necessary totake into account the effects of atom volumes and bondlengths, ρ redist has the unit of a charge density. Moredetails on its properties can be found in the SM [25]. FIG. 2. Charge density differences ρ diff of different mate-rials and their connection to adhesion values. Scales differfrom plot to plot. The effects of partial and full passivationof Fe(110) and Diamond(111) on ρ diff and W sep are shown. ρ redist values are shown as well and indicated by grey shadedareas. Units are the same as in Fig. 1. In the first row of Fig. 2, three materials with differentbonding types are depicted: van der Waals (vdW) bond-ing for double-layer graphene (a), abbreviated as Gr inall figures, metallic bonding for iron (b), and covalentbonding for diamond (c), abbreviated as C throughoutthe letter. The different scales in the color plots andthe line profiles, as well as the different values of E adh ,show a strong connection between E adh , ρ diff and in turn ρ redist . The adhesion increases by one order of magnitudebetween graphene and Fe, and between Fe and diamond,and a large increase in the magnitude of ρ diff can be ob-served as well. While the line profiles of ρ diff are similarin magnitude for iron and diamond, the metallic systemrearranges the charge much more uniformly at the in-terface center, while the charge is concentrated stronglyalong the directional carbon-carbon bonds in the insu-lator. While adhesion increases for C(111) compared toFe(110), ρ redist decreases. This indicates a different scal-ing of ρ redist for different bonding types, which we willrevisit in Fig. 3.Experiments have shown that lubricant additives con-taining sulfur or graphene can reduce the friction andwear of iron and steel significantly [23, 26, 27]. Like-wise, friction at diamond interfaces is greatly diminishedin the presence of H [28]. To investigate the root causeof these advantageous effects we investigate the influenceof these species on ρ redist and E adh at Fe and C inter-faces in the second row of Fig. 2. The leftmost panelresults from the addition of a / monolayer of sulfur atthe iron surfaces, which is the most favorable coverage forFe(110) [29]. This hinders charge accumulation at the in-terface compared to bare Fe, which in turn reduces E adh by a factor of ∼ .
5. The adhesion can be further reducedby higher coverage, as shown also in the case of phospho-rus [22]. In panel (e) we see that the iron surfaces arefully passivated by chemisorped layers of graphene, whichreduces the adhesion energy by one order of magnitudeinto the same range as double-layer graphene. Finally inpanel (f) we present results for diamond (111) surfaceswhere the presence of hydrogen termination leads to fullpassivation. In this case no covalent bonds are formedbetween the surfaces and ρ diff is reduced extremely atthe interface center. The adhesion for hydrogenated dia-mond is nearly two orders of magnitude smaller than forbare surfaces. Comparing (b) with (d) and (e), as well as(c) with (f), in Fig. 2, allows for a clear understanding ofthe lubrication mechanism of these passivating species,which consists of preventing charge accumulation at theinterface.In Fig. 3 we correlate ρ redist with the interfacial adhe-sion energy for a large set of different layered materials,metals and insulators/semiconductors. [30]It is remarkable that the redistribution of the chargeat the interface is so directly related to the adhesion ofsuch a wide variety of systems and surface terminations.The correlation within each bonding type is very good,with Pearson correlation coefficients for all fits are ∼ . FIG. 3. Adhesion energy versus charge redistribution ρ redist for (a) vdW bonded materials, (b) simple and noble metals(squares) and non-noble transition metals (circles), and (c)covalently bonded materials. C(Pan) specifies the Pandeyreconstruction. situation, something that usually is determined by anal-ysis of ∇ ρ and the Hessian matrix [31, 32]. For metals,Fig. 3(b), we see a distinct grouping of noble and simplemetals (squares) in the bottom left and the remainingtransition metals (circles) on the top right. This is ex-plained by the significant covalent bonding contributionof non-noble transition metals which increases both the E adh and ρ redist . It is important to note that for metalsand layered materials the linear relation holds also verywell if the independent variable ρ redist is replaced by theheight of the central peak of the planar average of ρ diff (see Fig. S1 in the SM [25]).In the following we analyze the potential corrugation∆ V . First we consider the relation between ∆ V andthe adhesion energy E adh . As can be seen in Fig. 4, wefind great correlation independent from the bond type.All data now gather around a single curve, which canbe very well fitted by a power law with exponent / ,∆ V = aE adh ( a =0.21 m J − ). Close relations betweenadhesion and friction have been discussed before [33–35],but a concrete functional relation has not been given yet.The added value of the power law is rather evident (al-though the exponent of / is not yet formally derived),as it permits to make precise predictions. The strong cor-relation between E adh and ∆ V is further evidence thatreducing adhesion by surface passivation using lubricantadditives leads to reduced friction (see Fig. 2), and isvery useful for the design of new lubricants, since E adh isusually more easily obtained by experiments than ∆ V .As a second step we analyze how ∆ V is related tothe charge redistribution occurring during sliding. Re-distribution of total charge [36] has been correlated be-fore to stacking fault energies of II-VI and III-V com-pounds (which are related to the PES) [37]. However,the partition of total charge into subsystems is some-what arbitrary [38] while our approach is based on the FIG. 4. PES corrugation ∆ V versus adhesion energy | E adh | .The main figure enlarges the region where most data are lo-cated while the inset shows all data. The black line fits alldata with a power law and their Pearson correlation coeffi-cient is 0 .
95. Symbols and colors are the same as in Fig. 3. charge density which is unambiguously defined and canbe evaluated using structure factors obtained from X-raydiffraction data with high accuracy [39, 40]. As a generaltrend we observe that for not ideal stacking the chargedensity is lower in the interface region than for the min-imum configuration. In Fig.5 we show this using chargedensity profiles for lateral configurations correspondingto the absolute minima and maxima of the PES for thesame materials as in Fig. 2(a)-(c) as well as Cu(001).As can be seen in the insets, the charge density profileat the center of the interface is lower for the maximumconfiguration than for the minimum configuration for allsystems.
FIG. 5. Total charge profiles for the interfaces of the materialspresented in the first row of Fig. 2 and Cu(001). Solid blacklines are for the minimum configurations; dashed red lines arefor the maxima. Charge is normalized to the respective valueat the center of the interface in the minimum configuration.
We can thus envisage an experiment where the simul-taneous measurement of ∆ V and ∆ ρ is performed andour finding is verified. Such experiment is schematicallyrepresented in Fig. 6. The interaction energy betweena probe moving on a crystalline substrate changes as a function of the relative lateral positions of the two bod-ies, as shown in Fig. 6(a) for copper on copper. We plot ρ in a slice at the center of the interface for three differentlateral configurations: the hollow [Fig. 6(b)], which is aminimum, an intermediate position [Fig. 6(c)], and theon-top configuration [Fig. 6(d)], which is a maximum.The arrow in Fig. 6(a) is visualizing the sliding path.We see that for the ideal fcc stacking the charge densityis high and is distributed rather uniformly [black solidcurve of Fig. 5(d)], while for the on-top configuration[Fig. 6(d)], significantly less charge [red dashed curve ofFig. 5(d)] is concentrated at the interface and it is lesshomogeneous than Fig. 6(b). This interfacial charge dif-ferences can be obtained, e.g., by a combination of Kelvinprobe force microscopy (KPFM), scanning tunneling mi-croscopy and atomic force microscopy (AFM) [41], or im-proved frequency modulated AFM [42]. KPFM is used tomeasure the work function φ locally [43, 44]. Changes of φ are directly related to charge density differences[45] atthe surface/interface [46]. It has also been shown that thework function of a material is correlated with its adhesivefriction [47]. For the example of the Cu(100) interface wecalculate a work function difference of ∆ φ = −
130 meVbetween two different surface stackings corresponding toFig. 6(b) and Fig. 6(d). The interaction potentials be-tween the tip and the surface, which for Cu(100) are rep-resented in Fig. 6(e), are also influenced by ∆ ρ and canbe measured by frequency shifts of an AFM [48]. FIG. 6. Potential energy surface of Cu(100) interface (a). To-tal charge densities at the interface for the minimum (b), aintermediate position (c), and the maximum (d), which aremarked on the sliding path shown by the arrow in (a). Inter-action potentials for minimum and maximum configurationsare shown in (e).
In summary, we have shown that the interfacial chargedensity and its variation during sliding are the basic phys-ical quantities, which determine adhesion, the PES, andthus the friction of a given atomically flat solid interface.The fact that the charge density is able to completelydefine the physical properties of a system is, of course, along-known result of density functional theory [49], buthere we have shown how to deduce important figures ofmerit on the tribological properties of an interface from ρ and ρ diff .The linear relationships between ρ redist , which quanti-fies the charge redistribution when an interface is formedfrom separated surfaces, and E adh is a very interestingresult that is relevant for the general field of interfacescience, beyond the tribological context in which it hasbeen presented here. We have shown that this simplelinear scaling holds for three different types of bonding:weak, physical vdW interactions, stronger metallic bond-ing with rather uniformly distributed charge, and verystrong directional covalent bonding. Simple power lawscaling of the PES corrugation ∆ V with E adh has alsobeen discovered where all investigated systems clusteraround a single curve.We have shown that the effectiveness of a certain classof friction reducing lubricant additives is to lower ρ redist ,which in turn leads to a reduced adhesion E adh and cor-rugation ∆ V .We have found a higher interfacial charge density ifthe investigated system is in an energy minimum thanin a maximum, and we connect this variation with thecorrugation of the PES. We have suggested that such aconnection can be experimentally observed by the simul-taneous measurement of electronic and frictional proper-ties during sliding.This work was partially supported by Materials De-sign at the Exascale (MaX) GA 676598 H2020 EINFRA-2015- 1, and the University of Modena and Reggio Emiliathrough the Fondi di Ateneo per la Ricerca (FAR) 2016project. M. C. Righi acknowledges support by the Eu-ropean Union through the MAX Centre of Excellence(Grant No. 676598). The authors thankfully acknowl-edge CINECA (Consorzio Interuniversitario del Nord estItaliano Per il Calcolo Automatico) for supercomput-ing resources through the project Italian SuperComput-ing Resource Allocation (ISCRA) B StressRx. Severalpictures in this Letter were created with the help ofXCrySDen [50]. ∗ [email protected][1] S. Y. Krylov and J. W. M. Frenken, physica status solidi(b) , 711 (2014).[2] T. D. B. Jacobs and A. Martini, Applied Mechanics Re-views , 060802 (2017).[3] Y. Zhou and J. Qu, ACS Applied Materials & Interfaces , 3209 (2017), pMID: 28029771.[4] H. Xiao and S. Liu, Materials & Design , 319 (2017).[5] I. Szlufarska, M. Chandross, and R. W. Carpick, Journalof Physics D: Applied Physics , 123001 (2008).[6] L. Prandtl, Zeitschrift f¨ur Angewandte Mathematik undMechanik , 85 (1928).[7] G. Tomlinson, Philosophical Magazine Series 7 , 905(1929).[8] W. Zhong and D. Tom´anek, Phys. Rev. Lett. , 3054 (1990).[9] D. Tom´anek, W. Zhong, and H. Thomas, Europhys.Lett. , 887 (1991).[10] M. C. Righi and M. Ferrario, Phys. Rev. Lett. , 176101(2007).[11] G. Zilibotti, M. C. Righi, and M. Ferrario, Phys. Rev.B , 075420 (2009).[12] G. Zilibotti and M. C. Righi, Langmuir , 6862 (2011),pMID: 21545120.[13] G. Levita, A. Cavaleiro, E. Molinari, T. Polcar, andM. C. Righi, The Journal of Physical Chemistry C ,13809 (2014).[14] S. Cahangirov, C. Ataca, M. Topsakal, H. Sahin, andS. Ciraci, Phys. Rev. Lett. , 126103 (2012).[15] M. Wolloch, G. Feldbauer, P. Mohn, J. Redinger, andA. Vernes, Phys. Rev. B , 195418 (2014).[16] M. Reguzzoni, A. Fasolino, E. Molinari, and M. C. Righi,Phys. Rev. B , 245434 (2012).[17] S. Cahangirov, S. Ciraci, and V. O. ¨Oz¸celik, Phys. Rev.B , 205428 (2013).[18] P. Restuccia and M. Righi, Carbon , 118 (2016).[19] L. Wang, X. Zhou, T. Ma, D. Liu, L. Gao, X. Li,J. Zhang, Y. Hu, H. Wang, Y. Dai, and J. Luo, Nanoscale , 10846 (2017).[20] A. R. Konicek, D. S. Grierson, P. U. P. A. Gilbert, W. G.Sawyer, A. V. Sumant, and R. W. Carpick, Phys. Rev.Lett. , 235502 (2008).[21] G. Zilibotti, S. Corni, and M. C. Righi, Phys. Rev. Lett. , 146101 (2013).[22] M. I. De Barros-Bouchet, M. C. Righi, D. Philippon,S. Mambingo-Doumbe, T. Le-Mogne, J. M. Martin, andA. Bouffet, RSC Adv. , 49270 (2015).[23] D. Marchetto, P. Restuccia, A. Ballestrazzi, M. Righi,A. Rota, and S. Valeri, Carbon , 375 (2017).[24] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcov-itch, Journal of Physics: Condensed Matter , 395502(2009).[25] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.121.026804 for:Computational details, additional plots of E adh versus ρ diff ( z = 0) and E adh versus 2 z × ρ diff , and ∆ V versus ∆ ρ . Furthermore detailes about the definition of ρ redist , a table with Pearson product-moment correlationcoefficients for all linear relations, and planar averagesof ρ diff for all investigated materials are given.[26] D. Berman, A. Erdemir, and A. V. Sumant, Carbon ,454 (2013).[27] M. C. Righi, S. Loehle, M. I. De Barros Bouchet,S. Mambingo-Doumbe, and J. M. Martin, RSC Adv. , 47753 (2016).[28] M.-I. De Barros Bouchet, G. Zilibotti, C. Matta, M. C.Righi, L. Vandenbulcke, B. Vacher, and J.-M. Martin,The Journal of Physical Chemistry C , 6966 (2012).[29] M. J. S. Spencer, I. K. Snook, and I. Yarovsky, The Jour-nal of Physical Chemistry B , 9604 (2005), pMID:16852156. [30] The planar averages of ρ diff of all investigated materialscan be found in the SM [25].[31] R. F. W. Bader, Atoms in Molecules - A Quantum Theory (University of Oxford Press, Oxford, 1990).[32] E. R. Johnson, S. Keinan, P. Mori-Snchez, J. Contreras-Garca, A. J. Cohen, and W. Yang, Journal of the Ameri-can Chemical Society , 6498 (2010), pMID: 20394428.[33] H. Yoshizawa, Y. L. Chen, and J. Israelachvili, The Jour-nal of Physical Chemistry , 4128 (1993).[34] K. Miyoshi, Tribology International , 605 (1999).[35] M. Urbakh, J. Klafter, D. Gourdon, and J. Israelachvili,Nature , 525 (2004).[36] J. C. Phillips and J. A. Van Vechten, Phys. Rev. Lett. , 1115 (1969).[37] S. Takeuchi, K. Suzuki, K. Maeda, and H. Iwanaga,Philosophical Magazine A , 171 (1985).[38] C. R. A. Catlow and A. M. Stoneham, Journal of PhysicsC: Solid State Physics , 4321 (1983).[39] T. S. Koritsanszky and P. Coppens, Chemical Reviews , 1583 (2001), pMID: 11709993.[40] H. Kasai, K. Tolborg, M. Sist, J. Zhang, V. R. Hathwar,M. Ø. Filsø, S. Cenedese, K. Sugimoto, J. Overgaard,E. Nishibori, and B. B. Iversen, Nature Materials ,249 (2018).[41] F. Mohn, L. Gross, N. Moll, and G. Meyer, Nature Nan- otechnology , 227 (2012).[42] F. Albrecht, J. Repp, M. Fleischmann, M. Scheer,M. Ondr´aˇcek, and P. Jel´ınek, Phys. Rev. Lett. ,076101 (2015).[43] S. Sadewasser, T. Glatzel, M. Rusu, A. J¨ager-Waldau,and M. C. Lux-Steiner, Applied Physics Letters , 2979(2002).[44] T. Filleter, K. V. Emtsev, T. Seyller, and R. Bennewitz,Applied Physics Letters , 133117 (2008).[45] The connection is: ∆ φ = − e(cid:15) ∆ p , whith the changein the surface dipole ∆ p depending on ρ diff , ∆ p ( z ) = (cid:82) zρ diff ( z ) dz [51].[46] S. Ling, M. B. Watkins, and A. L. Shluger, Phys. Chem.Chem. Phys. , 19615 (2013).[47] Y. Li and D. Y. Li, Journal of Applied Physics , 7961(2004).[48] N. M. Caffrey, K. Buchmann, N. Hauptmann, C. Lazo,P. Ferriani, S. Heinze, and R. Berndt, Nano Letters ,5156 (2015), pMID: 26218345.[49] P. Hohenberg and W. Kohn, Physical Review , 864(1964).[50] A. Kokalj, Computational Materials Science , 155(2003), proceedings of the Symposium on Software De-velopment for Process and Materials Design.[51] T. C. Leung, C. L. Kao, W. S. Su, Y. J. Feng, and C. T.Chan, Phys. Rev. B68