Modulating Bond Lengths via Backdonation: A First ... - ACS Publications

Mar 4, 2016 - Department of Physical and Life Sciences, Alfred State College, SUNY ... Department of Chemistry, State University of New York at Buffal...
0 downloads 0 Views 3MB Size
Article pubs.acs.org/JPCC

Modulating Bond Lengths via Backdonation: A First-Principles Investigation of a Quinonoid Zwitterion Adsorbed to Coinage Metal Surfaces Scott Simpson,† James Hooper,*,‡ Daniel P. Miller,§ Donna A. Kunkel,∥ Axel Enders,∥ and Eva Zurek*,§ †

Department States ‡ Department § Department ∥ Department States

of Physical and Life Sciences, Alfred State College, SUNY College of Technology, Alfred, New York 14806, United of Theoretical Chemistry, Faculty of Chemistry, Jagiellonian University, R. Ingardena 3, 30-060 Krakow, Poland of Chemistry, State University of New York at Buffalo, Buffalo, New York 14260-3000, United States of Physics and Astronomy, University of NebraskaLincoln, 855 N 16th Street, Lincoln, Nebraska 68588, United

S Supporting Information *

ABSTRACT: First-principles calculations reveal that upon adsorption to the Cu(111) surface, the C−C single bonds within the p-benzoquinonemonoimine zwitterion (ZI) contract by about 6%. A detailed analysis reveals that the bond shortening is primarily a result of backdonation from Cu orbitals of s and d symmetry to the lowest unoccupied orbital (LUMO) of the ZI. This LUMO is π*-antibonding across the molecule and π-bonding across the C−C bond that shortens. We illustrate that the level alignment between the Fermi level of the surface and the frontier molecular orbitals of the ZI, the topology of the LUMO, and the distance between the substrate and the adsorbate are important factors enabling bond strengthening via backdonation. An extended transition state−natural orbitals for chemical valence (ETS-NOCV) analysis is applied to molecular models for this system, and it confirms that the surface → LUMO backdonation on Cu(111) is larger than on Ag(111) and Au(111).



INTRODUCTION Designing organic networks whose molecular architecture gives rise to specific properties requires meticulous fabrication by precisely manipulating the molecules as well as the interactions between them. The reactivity, aromaticity, and gap between the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) of individual molecules, and of the self-assembled networks they form, are intimately linked to their geometries, which can be modulated by the surface via the electronic interactions that occur at the metal− organic interface.1 For example, a growing number of novel synthesis strategies that exploit on-surface thermal activation have been developed.2−8 The changes in the bond lengths that occur when carbon monoxide adsorbs to the metal surface can be explained in terms of Blyholder’s model.9 Because of the electron donation from the HOMO of CO to the d-bands on the metal surface with a concomitant backdonation from the d-bands to the LUMO, the C−O bond stretches. In many cases, the electronic reorganization at the substrate−adsorbate interface can be understood via the Blyholder picture adapted to unsaturated hydrocarbons: the Dewar−Chatt−Duncanson (DCD) model.10,11 The DCD model has been applied to various © 2016 American Chemical Society

organics adsorbed on metallic substrates including ethylene, benzene,12,13 acetylene,14 porphine,15 tetraphenylporphyrin,16 and the quinonoid zwitterion (ZI) p-benzoquinonemonoimine illustrated in Figure 1a.17,18 In both the adsorption of organic molecules on surfaces and the interaction of metal atoms with ligands, the electron density is oftentimes donated into an orbital that is antibonding.19 Because of this, the metal−organic interaction results in a lengthening of the bond into which electron density is donated. For carbon monoxide, stronger substrate−adsorbate interactions yield a longer C−O bond.20 In some circumstances, the nodal structure of the LUMO of the adsorbate or ligand may be bonding across one of the bonds. Consider butadiene as an example, wherein the HOMO is π-bonding across C1−C2/ C3−C4 and π*-antibonding across C2−C3, and the LUMO is π-bonding across C2−C3 and π*-antibonding across C1−C2/ C3−C4 (Figure 1b). When this ligand coordinates to a metal center in its cisoid conformation, the ensuing charge transfer can yield complexes in which the C2−C3 bond is shorter than Received: January 12, 2016 Revised: March 1, 2016 Published: March 4, 2016 6633

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C

geometric changes, nor discussed how they are linked to the electronic perturbations at the metal−organic interface. We find that the length of the 2 C−C bond in the ZI from Figure 1a can be tuned by the choice of the substrate. Whereas a pronounced shortening of about 6% (0.098 Å) occurs when the molecule is adsorbed on Cu(111), on Au(111) the contraction is only 0.016 Å. Energy decomposition schemes such as the ETSNOCV analysis used by Braunschweig are commonly used for molecules, but their application to molecular models for surface-adsorbed systems is complicated by the large number of metal orbitals that can participate in the electron transfer between the substrate and the adsorbate. Herein, we apply the ETS-NOCV method for the first time to such a system, thereby illustrating that the bond length changes can be understood in the degree of forward (molecule → surface) and back (surface → molecule) donation. The molecular orbital into which electron density is donated in the ZI is π*-antibonding across the molecule and π-bonding across the C−C bond that shortens, and calculations on a molecular model system suggest that the Cu orbitals involved in the backdonation are of s and d symmetry.

Figure 1. (a) A ball and stick and skeletal drawing of the parent quinonoid ZI (p-benzoquinonemonoimine compound from the class of N-alkyldiaminoresorcinones). The HOMO (red/blue) and LUMO (purple/green) of the optimized ZI. (b) A ball and stick and skeletal drawing of cis-1,3-butadiene. The isovalue plots of the HOMO and the LUMO of butadiene. (c) A ball and stick and skeletal drawing of the [Pt(PMe3)2(B2Ph2)] compound studied in ref 25. The B2Ph2 ligand and its LUMO. The orbital isosurfaces are illustrated using an isovalue of ±0.03 au (a, c) and ±0.04 au (b). (d) Lewis structures of the neutral ZI and the ZI with a charge of −2.



RESULTS AND DISCUSSION Periodic DFT calculations that accounted for dispersion semiempirically (PBE-D2, revPBE-D3) and those that approximate dispersion using a nonlocal correlation functional (vdWDF2, vdW-DF-optB88) were used to optimize models for the isolated ZI shown in Figure 1a on Cu(111). The interaction energies (ΔE) were calculated as being −0.73, −1.80, −2.17, and −2.37 eV via the vdW-DF2, vdW-DF-optB88, PBE-D2, and revPBE-D3 functionals, respectively. These trends are expected since vdW-DF2 typically underbinds and the Grimme corrections often overbind adsorbate molecules to metal surfaces, while vdW-DF-optB88 tends to be intermediate, often yielding results that compare better with experiment.30−32 Because the computational cost of the vdW-DF-optB88 functional depends on the number of FFT grid points, it becomes expensive for large unit cells, whereas the cost of revPBE-D3 and PBE-D2 is significantly lower. Because Grimme’s D2 correction does not have parameters for gold, and considering the reasonable agreement between the binding energies and bond length changes (see Figure 2) computed with vdW-DF-optB88 and revPBE-D3, the latter functional was employed for the calculations presented in this article, unless otherwise stated. The ΔEs computed for the adsorption of the ZI on the Cu(111), Ag(111), and Au(111) surfaces were −2.37, −1.60, and −1.83 eV, respectively. Because of the tendency of bonds to elongate upon adsorption, we were quite surprised when comparing the gas phase geometric parameters of the ZI with those of the molecule optimized on Cu(111). Figure 2 shows that the C−O, C−N, 1 C−C, and 3 C−C bonds elongated upon surface adsorption, as expected (see Figure 1a for the bond labels). However, the 2 C−C bond shortened significantly for all of the functionals employed. Upon further analysis, we realized that this behavior can be rationalized. As illustrated in Figure 1a, the LUMO of the zwitterion possesses a node that is coincident with the mirror plane passing through the center of the molecule and is of C−O and C−N π*-character. However, this LUMO is also π-bonding across the 2 C−C bond. Therefore, backdonation into this orbital would explain the bond contraction. The 2 C−C bond shortened on all three of the coinage metals studied, but the effect became less pronounced

that of the C1−C2/C3−C4 bonds, or to situations where the three sets of bonds are roughly equal in length.21,22 Adsorption of cis-1,3-butadiene to metal surfaces can result in geometries where only the terminal C−C bonds have been elongated or ones where the C2−C3 bond is also somewhat stretched.23,24 In a recent study on the discrete transition-metal π-diborene complex illustrated in Figure 1c, Braunschweig’s group demonstrated, for the first time, that π-backdonation can lead to bond strengthening in a transition metal π-diborene complex.25−27 Density functional theory (DFT) calculations showed that an important component of the bonding between [(Me3P)2M] and B2Ph2 (Me = methyl, Ph = phenyl, and M = Pd, Pt) is M → B2 backdonation. The empty π-acceptor orbital was localized across the B−B bond of B2Ph2, and it was oriented orthogonal to the MP2 plane (π⊥ symmetry). An extended transition state−natural orbitals for chemical valence (ETS-NOCV) analysis illustrated that the π-backbonding interaction was strong and stabilizing. The B−B bond lengthened by only 0.9−2% as a result of coordination to the metal center, much less than the C−C lengthening of 7.7% in comparable systems. This theoretical prediction inspired experimental work that led to the synthesis of an analogous complex with the same type of bonding, [Pt(PEt3)2(B2(Dur)2)] (Dur = duryl = 2,3,5,6-tetramethylphenyl), and has sparked the computational search for other transition-metal complexes28 and silylenes29 that can be employed to shorten bonds in diborane and π-diborenes. Herein, we present the results of detailed DFT investigations of the interaction of the ZI illustrated in Figure 1a,d with the coinage metal surfaces Cu(111), Ag(111), and Au(111). The nodal structure of the LUMO is similar to that of butadiene in the sense that it is π-bonding across adjacent carbon atoms. We have previously illustrated that the bonding of these surfacesupported systems can essentially be explained by the DCD model17,18 but have not commented on the intramolecular 6634

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C

copper, 0.95 eV on silver, and 1.10 eV on gold. The largest reduction in the gap occurred for the most distorted geometry, as expected. When the ZI and the surface slabs are brought together, their wave functions overlap and the states containing character from the HOMO and LUMO of the ZI broaden (see Figure 3). Nevertheless, the energy differences between the middle of the HOMO-like and LUMO-like states correlate well with the gaps of the isolated, distorted ZIs. The difference in the magnitude of the 2 C−C bond contraction on the three coinage metals can be explained by considering where the ZI’s LUMO lies with respect to EF. Because the Fermi level of copper is about 1.5 eV higher than that of silver and gold, on Cu(111) a majority of the LUMO-like states are found below EF, whereas on Ag(111) the LUMO straddles EF, and on Au(111) the center of the peak in the DOS exhibiting character from the LUMO lies ∼0.5 eV above EF. All of the systems exhibit HOMO → surface forward donation, but the amount of surface → LUMO backdonation decreases significantly in going from copper to gold. Thus, the PDOS confirm that the interaction of the ZI with the (111) surfaces of the coinage metals can, for the most part, be explained by the DCD model, and backdonation from the metal into the LUMO of the adsorbate, which is π-bonding across the 2 C−C bond, leads to a shortening of this bond. The bond contraction correlates with the amount of backdonation, which is largest on Cu(111) and smallest on Au(111). Our calculations have suggested that the positions of the atoms comprising strongly activating and deactivating functional groups in the adsorbate, with respect to the underlying surface atoms, are one of the key factors determining the most stable adsorption sites for benzene derivatives on Ag(111).33 The most stable binding site we found for ZI adsorption to Cu(111) placed the center of the molecule directly over a Hhcp site, with the nitrogen atoms being close to a top (T) site and the oxygen atoms over Hfcc sites.34 This configuration is illustrated in the inset in Figure 4b. The notation we use to describe it, Hhcp0, provides the location of the center of the carbon ring and the number of degrees it has been rotated relative to the ⟨112⟩ direction. The interaction energies, geometric distortions, and PDOS plots discussed above were computed for the most stable configurations. Could it be that the position of the ZI on the Cu(111) surface affects the 2 C−C contraction? In order to test this, we optimized systems where the location of the zwitterion was varied (see the Supporting Information for more details), which resulted in a computed corrugation energy of 0.54 eV. Figure 4a provides the bond length changes for the most stable (Hhcp0) and least stable (T0) configurations we found as well as one that was intermediate in energy (T30). The relative interaction energies were calculated as being 0/0.11/0.54 eV, and the 2 C−C bond length decreased by 0.096, 0.098, and 0.066 Å on the Hhcp0, T30, and T0 sites, respectively. Our results show that larger binding energies generally result in a greater contraction of the 2 C−C bond. On Ag(111) and Au(111), the corrugations in the potential energy surface were calculated as being at least 0.17 and 0.15 eV, respectively, and the 2 C−C bonds shrunk by 0.050 and 0.016 Å, respectively, at the most stable adsorption sites. The Supporting Information contains more details about the specific binding sites tested for Cu, Ag, and Au. The ZIs self-assemble to form zipper-like rows on the Cu(111) surface, as shown in our calculated structure model (Figure 4b) and the STM image in Figure 4c. Because the

Figure 2. Deviation of bond lengths of the zwitterion adsorbed on Cu(111) compared to the gas phase calculated with different functionals that account for dispersion, as reported in Å. A negative/ positive value denotes a shortening/elongation of the bond upon adsorption. The bond labels are defined in Figure 1. The vdW-DF2, vdW-DF-optB88, revPBE-D3, and PBE-D2 calculations were performed using a coverage of 0.612 molecules/nm2.

for the heavier metals of Ag and Au. The contraction was calculated as being 0.098, 0.050, and 0.016 Å on the (111) surfaces of Cu, Ag, and Au, respectively. To explain this trend, we calculated the projected density of states (PDOS) (see Figure 3) and analyzed the alignment of the adsorbate’s frontier orbitals with the Fermi levels (EF) of the metallic substrates. Even though the revPBE-D3 nonhybrid functional employed in this study underestimates band gaps, it is instructive to examine how the ZI’s gas phase HOMO−LUMO gap of 1.3 eV is affected by the geometric distortion it experiences upon adsorption. Single-point calculations on isolated molecules in the geometry adopted on the surface yielded gaps of 0.40 eV on

Figure 3. PDOS of the isolated ZI and of the ZI adsorbed to the Cu(111), Ag(111), and Au(111) surfaces with revPBE-D3 optimized lattice constants. The total PDOS of the ZI is indicated by the black line. The blue line corresponds to the PDOS of the H and C atoms that belong to the σν symmetry plane. These atoms contribute toward the HOMO of the ZI, but not to the LUMO (see Figure 1a). The red line corresponds to the PDOS of the carbon atoms that are directly bonded to the amino groups. These atoms contribute toward the LUMO of the ZI, but not to the HOMO. In all plots EF has been adjusted to zero. The coverage was 0.613/0.466/0.456 molecules/nm2 for Cu/Ag/Au, respectively. 6635

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C

calculations were carried out to decompose the BSSE (basis set superposition error) uncorrected interaction energies into chemically intuitive terms such as the energy associated with distorting the ZI into the geometry it assumes postadsorption (ΔEgeo), the repulsive steric energy (ΔEsteric), which includes the Pauli repulsion and the electrostatic interactions between the interacting fragments, the attractive interactions due to dispersion (ΔEdisp), and the orbital interactions (ΔEoi), which relates to the charge transfer between fragments.36 For the Cu(111) surface, we considered the least stable, T0, and intermediate, T30, sites for comparison with the most stable, Hhcp0, site. The molecular revPBE-D3 binding energies were within 0.1−0.3 eV of the results obtained with the revPBE-D3 periodic calculations. The small differences can be attributed in part to the different basis sets (atom-centered vs plane-wave) and the finite surface models. The BSSE corrections to the binding energies were found to be less than 13% of the total interaction energy (0.29, 0.14, and 0.13 for Cu, Ag, and Au, respectively; see the Supporting Information). Figure 5a,b compares the non-BSSE corrected binding energies, and the bars are colored to show the relative weights of the components that comprise this binding energy. It is instructive to correlate the trends in the contributions to the total binding energy for the three different sites on Cu(111) with their geometries and the change in the occupation of the HOMO and the LUMO upon adsorption (see Figure 5a(i−ii)). Shorter surface−adsorbate distances resulted in a larger magnitude of ΔEsteric. ΔEdisp was found to be relatively independent of the adsorption site. The closer proximity of the substrate and adsorbate led to a larger overlap of their wave functions and a more stabilizing ΔEoi. In general, we found that a larger magnitude of the orbital interaction energy gave rise to a shorter 2 C−C bond length over a wide range of cluster models (see the Supporting Information). The reason for this is that the increased orbital interactions led to a significantly larger surface → LUMO backdonation, as illustrated by the estimated orbital occupations in Figure 5a(i−ii). This is in line with the proposition of Michaelides et al. that the preferred binding sites are driven by a tendency to maximize the overlap between the surface and adsorbate wave functions.37 The HOMO → surface electron transfer, on the other hand, did not differ much for the three adsorption sites. A comparison of the terms comprising the binding energies of the ZI adsorbed to the most stable binding site on Cu(111), Ag(111), and Au(111) in Figure 5b(i) further illustrates that the absolute magnitude of the terms comprising the binding energies are the largest on Cu(111) because it has the shortest substrate−adsorbate distance and the largest geometric distortion of the ZI. The estimated occupations of the HOMO and LUMO in Figure 5b(ii) confirm that the degree of surface → LUMO backdonation is significantly larger on Cu(111) than on the heavier metals, and this is the reason for the much shorter 2 C−C distance. The depopulation of the HOMO upon adsorption on the three different metals is not as sensitive to the identity of the metal as is the population of the LUMO. To help better quantify the contributions to ΔEoi on the three substrates and how they relate to the ZI molecular orbitals, we employed the ETS-NOCV energy decomposition scheme.38 ETS-NOCV was used by Braunschweig to show that π-backbonding stabilizes the [Pt(PMe3)2(B2Ph2)] complex by 4.16 eV.25 It has been applied to unravel the bonding in a wide range of systems including hydrogen bonding39 in guanine

Figure 4. (a) As in Figure 2, except for the revPBE-D3 functional with the center of the ring comprising the ZI positioned on different sites of the Cu(111) surface. The changes in the bond lengths on the Hhcp0, T30, and T0 sites (the Hhcp0 is illustrated in the inset of (b)) were computed with a coverage of 0.612 molecules/nm2. The ZI−surface distance (calculated by subtracting the z-coordinate of the atoms in the top of the unoptimized surface and the average of the z-coordinate of the C, O, and N atoms in the ZI optimized on the surface) are 2.29 Å (Hhcp0), 2.38 Å (T30), and 2.53 Å (T0). The high-coverage (HC) data were obtained for the network in (b) that was constructed as a model for the STM image of the ZI on Cu(111) in (c). The HC network had a coverage of 1.36 molecules/nm2 and a molecular spacing of 1.25 molecules/nm in the 2D chains.

surface coverage can influence the strength of the substrate− adsorbate interaction, and hence the charge transfer at the metal−organic interface,18,35 we investigated its effect on the 2 C−C shortening. The best model we previously found for the zipper-like rows, shown in Figure 4b, placed each molecule at either the Hhcp0 site or the analogous Hfcc180 site, and the spacing between molecules, which is a function of the coverage,34 was 1.25 molecules/nm. This model agrees well with the experimentally observed ZI chains as it reproduces the intermolecular spacing and the alignment of the chains along three distinct crystallographic directions, which are the Cu ⟨110⟩ directions. Because of the strong interaction between the ZI and Cu(111), the bond length changes observed in this high-coverage model agreed well with those computed at lower coverages (see Figure 4a), where the model for the zipper-like rows is labeled HC. On the other hand, our previous calculations18 suggest that strong hydrogen bonding between the ZIs in the networks formed on Ag(111) and Au(111) results in a weakening of the substrate−adsorbate interactions and a decrease in the charge transfer. Since the 2 C−C bond contraction on these metals at low coverages was not nearly as pronounced as on Cu(111), it may be that at the coverages observed in experiment there is very little 2 C−C shortening, especially on Au(111). Single point molecular calculations on the most stable ZI/ surface configurations found in our revPBE-D3 periodic 6636

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C

number of NOCV pair contributions with small, near negligible, energies. However, the visualization of the NOCV pairs with the largest contributions, and how much each pair contributes to the total ΔEoi, is instructive. This is what we show in Figure 5c, wherein each colored block represents a single NOCV pair. The color of the block describes the direction of charge flow. The intent behind using ETS-NOCV here is to further quantify the amount of HOMO → surface and surface → LUMO donation on the three metal clusters. First, we consider the ZI−Cu(111) interaction. The flow of charge that is described by the first and most significant NOCV, labeled 1 in Figure 5c, corresponds very well, both qualitatively and quantitatively, with charge transfer from the surface into the ZI LUMO. The interaction energy associated with this “first” LUMO-based pair is −2.54 eV, which accounts for ∼31% of ΔEoi. The fragment orbital analysis suggests that the population of the ZI LUMO associated with this NOCV is 1.44e, which, yet again, is consistent with the large occupancy of the LUMO in Figure 5b(ii). The second highest energy NOCV pair corresponds mostly to a charge depletion of the ZI HOMO (∼0.15e) and contributes −0.56 eV to ΔEoi. The rest of the NOCVs that are colored for Cu in Figure 5c correspond to charge transfer from ZI to the surface cluster, but they all involve molecular orbitals which are lower in energy than the HOMO; for example, the third NOCV involves both the HOMO and the HOMO−10. To simplify the comparison with Ag and Au, however, let us focus on how the LUMO and HOMO orbitals contribute to ΔEoi. In total, 0.29e of charge density is associated with HOMO → Cu(111) donation, and 1.48e is involved in Cu(111) → LUMO backdonation. The highest energy NOCV pair from our model Ag(111) system also corresponds to charge donation into the ZI LUMO and contributes 0.41 eV (roughly 20%) to ΔEoi, which is substantially less than on Cu(111). The highest energy NOCV pair for Au(111), on the other hand, corresponds to HOMO → Au(111) charge transfer; it contributes 0.71 eV (roughly 33%) to ΔEoi. The third and fourth most significant NOCVs from the Au(111) analysis correspond to Au(111) → LUMO charge transfer, and they contribute 0.25 eV to ΔEoi, i.e., 12% of ΔEoi. To summarize, the parts of ΔEoi that are associated with backdonation from the surface into the LUMO were calculated to be −2.54, −0.41, and −0.25 eV for the copper, silver, and gold models, respectively. This agrees quite well with the estimated occupation of the ZI LUMO in each system: 1.48e on Cu(111), 0.56e on Ag(111), and 0.07e on Au(111). It also shows that there is less charge transferred into the LUMO on Ag(111) and Au(111) and accounts for the largest single difference in the NOCV contributions to ΔEoi. Therefore, the pronounced contraction of the 2 C−C bond on Cu(111) and its smaller contraction on the other metals can be attributed to the amount of surface → LUMO charge transfer. Thus, the NOCV analysis confirms the conclusions we reached based upon analyzing the PDOS plots in Figure 3. The C−C bond lengths within ethane, benzene, and ethene were calculated as being 1.538, 1.404, and 1.339 Å. The length of the 2 C−C bond of the ZI in our molecular calculations (1.552 Å) is typical of a single bond, but all of the carbon− carbon bonds in the ZI adsorbed on Cu(111) are intermediate between a single and a double bond. This made us wonder if bond contraction resulting from surface → LUMO backdonation could impact the zwitterion’s propensity for

Figure 5. (a) (i) A bar graph comparing the decomposition of the BSSE uncorrected interaction energy, ΔEbind, of ZI at the Hhcp0, T30, and T0 binding sites on Cu(111). ΔEbind = ΔEgeo + ΔEsteric + ΔEoi + ΔEdisp, and the colors correspond to the relative weight that each term contributes to ΔEbind. All values are given in eV, and 114-atom clusters were used to model the surface in the molecular revPBE-D3/TZ2P/ TZ2P+ calculations. (ii) A bar graph comparing the gross populations of the LUMO and HOMO of the ZI when adsorbed to the Hhcp0, T30, and T0 binding sites. (b) (i) + (ii) Same as in part (a), but for the most stable binding sites found for Cu, Ag, and Au. The ZI−surface distances were calculated as being 2.29, 2.75, and 2.84 Å on Cu(111), Ag(111), and Au(111). (c) A schematic of how ETS-NOCV decomposes the ΔEoi term for the ZI adsorbed to the Cu, Ag, and Au models. Each colored block corresponds to the energy contribution of an NOCV (the eight highest in energy are shown), and the color corresponds to the type of charge transfer event that best describes the NOCV. Blue: surface→ LUMO; red: involves various occupied ZI MOs including the HOMO; yellow: has features of both blue and red.

quartets40 as well as the adsorption of formaldehyde,41 ethene,42 and ethylene42,43 to Cu(I) and Ag(I) sites. This scheme decomposes ΔEoi into contributions from complementary orbital pairs (NOCVs) wherein each NOCV pair defines a distinct charge transfer event between the substrate and the ZI in the geometry it assumes on the surface.38 The ETS-NOCV scheme for extended systems is currently being developed. This article presents the first application of the ETS-NOCV method to an organic molecule adsorbed to a metallic cluster. We classify the NOCV pairs, qualitatively, by looking at the flow of charge density associated with each pair and, quantitatively, by using the change in the ZI fragment orbital populations that are printed as part of the fragment orbital analysis of each NOCV that is implemented in ADF. We note that the nature of the finite cluster tends to result in a large 6637

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C

Figure 6. (a) An approximate interaction diagram for the formation of the molecular model composed of the ZI and four copper atoms (ZI−4Cu). (b) Isosurfaces, with isovalues ranging from 0.02 to 0.03 au, of MOs containing character from the LUMO of the ZI.

aromaticity? To answer this question, we calculated the zz component of the nucleus independent chemical shift44,45 tensor 1 Å above the center of the molecular ring (NICS(1)zz), as this has been found to be an excellent indicator of πaromaticity in planar molecules.46 The NICS(1)zz computed at the revPBE-D3/TZP level of theory for a gas phase benzene molecule (−27.8 ppm) and a free ZI molecule (+11.5 ppm) are typical for aromatic and antiaromatic systems, respectively. We utilized a simple molecular model for the adsorbed ZI that has two Cu atoms interacting with the nitrogens within the amine groups and two more Cu atoms interacting with the oxygen atoms on the carbonyl groups (we refer to this system as ZI− 4Cu). The 1 C−C, 2 C−C, and 3 C−C bond lengths in the optimized system were 1.406, 1.440, and 1.397 Å, respectively, in close agreement with our periodic calculations. The NICS(1)zz value for this model system was calculated to be −12.4 ppm, showing that the interaction with the Cu atoms transmogrifies this molecule from an antiaromatic into an aromatic system. The NICS values we calculated at the center of the rings were −7.5, + 7.1, −9.1, and −8.2 ppm for benzene, the ZI, the ZI dianion, and the ZI interacting with four copper atoms, respectively. This computational experiment suggests that it may be possible to find molecular transition metal containing complexes that can interact with the ZI and induce bond-length changes similar to those predicted to occur on the Cu(111) surface. An additional computational experiment where we optimized the geometry of the same four metal atom cluster interacting with the ZI, but replaced the copper atoms with platinum atoms (we refer to this system ZI−4Pt), resulted in 1 C−C, 2 C−C, and 3 C−C bond lengths of 1.439, 1.402, and 1.393 Å. The 2 C−C bond in this model system is 0.038 Å shorter for ZI−4Pt than it is for ZI−4Cu, suggesting that the interaction of the ZI with either a platinum surface or a platinum containing molecular complex may also lead to this type of bond shortening. Finally, we calculated an approximate MO interaction diagram for the formation of the ZI−4Cu system from four noninteracting Cu atoms and the ZI, shown in Figure 6a.

Isosurfaces of the occupied MOs of this organometallic complex that exhibit character arising from the LUMO of the ZI, the 12a″, 13a″, 14a″, 15a″, and 22a″, are illustrated in Figure 6b. The first five of these MOs are composed primarily of Cu dstates (82.0−94.6%) and contain only between 1.2 and 11.6% character from the LUMO of the ZI. The 22a″ MO, on the other hand, is composed of the ZI LUMO (38.8%), Cu s-states (36.4%), and Cu d-states (14.3%). The 23a″ LUMO contains character arising from the ZI LUMO as well as Cu s, p, and d states. Thus, in this model system, the backdonation into the ZI-LUMO occurs not only from the Cu d but also from the Cu s states. We note that previous calculations have found differences between Blyholder’s model of CO adsorption to the (111) surfaces of Pt, Cu, and Al, in particular showing that the surface s and p bands are of prime importance for this interaction.47 Figure 1c provides a schematic Lewis dot structure of the ZI assuming a full transfer of two electrons from Cu into its LUMO. The bond strengthening via backdonation on a metal surface is not unique to this quinonoid zwitterion. Indeed, it is found for the adsorption of cis-1,3-butadiene on certain surfaces, and we predict a similar bond shortening for another organic compound, benzoquinone. Our preliminary investigations suggest the regiochemistry (ortho-, meta-, and para-) of the compound affects the importance of this bond shortening mechanism (see Supporting Information).



CONCLUSIONS Herein, we illustrate how the choice of the substrate may be used to modulate the degree of bond strengthening via πbackbonding in an organic adsorbate, thereby influencing the geometric parameters and properties of the adsorbed molecule. Our density functional theory (DFT) calculations show that backdonation into the LUMO of the p-benzoquinonemonoimine zwitterion (ZI) shortens the single C−C bond, so that all of the C−C bonds in the molecule approach the length of those found within benzene. The effect is pronounced on the Cu(111) surface because most of the states arising from 6638

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C

frozen. Scalar relativistic effects were included by using the zeroth-order regular approximation (ZORA) Hamiltonian60 and ZORA basis sets. The Cu(111), Ag(111), and Au(111) surfaces were modeled using a 114-atom cluster consisting of three layers with 48/38/38 atoms in the top/middle/bottom layers. The geometries of the slabs were obtained from our periodic VASP calculations. The basis set superposition error (BSSE) was obtained using the counterpoise method.

the LUMO of the molecule fall below the Fermi level of the substrate. The degree of bond contraction decreases in going to Ag(111) and Au(111) because the degree of backdonation decreases. The LUMO of the ZI is π*-antibonding across the molecule but π-bonding across the C−C bond that shortens, and other molecules, such as benzoquinone, are predicted to behave similarly. An ETS-NOCV (extended transition state− natural orbitals for chemical valence) analysis was used to meticulously analyze the different contributions to the orbital interaction energy in molecular models for the ZI/surface systems, confirming that the degree of bond contraction is related to the degree of surface → LUMO backdonation. This article represents the first application of the ETS-NOCV analysis to cluster models for molecules adsorbed to surfaces.



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.6b00360. Optimized coordinates for each system, selected interatomic distances, an in-depth ETS-NOCV analysis for select systems, in addition to an analysis of the binding sites (PDF)



COMPUTATIONAL METHODS Periodic DFT calculations (geometry optimizations, electronic densities of states, and charge densities) were performed using the Vienna Ab-initio Simulation Package (VASP) version 5.3.5.48 The projector augmented wave (PAW) method49 was used to treat the core states along with a plane-wave energy cutoff of 500 eV, and the C/N/O 2s/2p, H 1s, Cu 4s/3d, Ag 5s/ 4d, and Au 6s/5d electrons were treated explicitly. The following functionals were employed to account for dispersion interactions:32 the GGA functional of Perdew−Burke− Ernzerhof50 with the DFT-D2 method of Grimme51 (PBED2), the revised PBE functional with Grimme’s zero-damping DFT-D3 method52 (revPBE-D3), the optimized Becke8853 van der Waals functional (vdW-DF-optB88) developed by Klimeš et al.,54−57 and the nonlocal correlation functional of Langreth and Lundqvist (vdW-DF2).54,58 The Γ-centered Monkhorst− Pack scheme was used to generate k-point grid meshes, and we employed a 7 × 7 × 1 k-mesh, which gave rise to energies that were converged to within 2 meV/atom of energies computed with a larger k-mesh. The geometries were converged such that the magnitude of the largest force acting on the optimized atoms was less than 0.05 eV/Å. The “free” zwitterion was optimized in a box measuring 20 × 20 × 20 Å. DOS plots were obtained using a 15 × 15 × 1 k-mesh. A dipole correction was applied along the direction perpendicular to the metal surface using the LDIPOL tag as implemented in VASP. Surfaces models composed of four layers were employed, and a vacuum spacing of ∼30 Å was used to separate two slabs. In the geometry optimizations, the molecular coordinates and upper two layers of the surface were allowed to freely relax, whereas the bottom two layers were kept fixed. For the calculations of the ZI/Cu(111) system computed with different functionals, the experimental bulk lattice constant of 3.615 Å was used for Cu, and the surface was simulated using a 100 atom slab. The remainder of the computations were carried out with the optimized revPBE-D3 lattice constants of 3.542, 4.103, and 4.138 Å for Cu, Ag, and Au, respectively, and a 120-atom simulation cell. In the Supporting Information we show that the revPBE-D3 interaction energies and bond lengths of the ZI were relatively insensitive to the lattice constant employed (theoretical vs experimental) for the Cu surface. Molecular calculations were carried out using the Amsterdam density functional (ADF) software package36,59 and the revPBE-D3 functional. The basis functions on all of the atoms consisted of a triple-ζ Slater-type basis set with polarization functions (TZ2P+ for Cu, TZ2P for everything else) from the ADF basis set library. The core−shells up to 1s for C, N, O, 3p for Cu, 4p for Ag, and 4f for Au were kept



AUTHOR INFORMATION

Corresponding Authors

*E-mail [email protected] (J.H.). *E-mail ezurek@buffalo.edu (E.Z.). Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We acknowledge the National Science Foundation (NSF) (No. HRD-1435163) for financial support. E.Z. thanks the Alfred P. Sloan Foundation for a Research Fellowship (2013-2015). J.H. acknowledges financial support from the Homing Plus Program (HOMING PLUS/2012-6/4) granted by the Foundation for the Polish Science and cofinanced by the European Regional Development Fund. A.E. acknowledges support from the National Science Foundation through the Materials Research Science and Engineering Center (Grant DMR-1420645). Support from the Center of Computational Research at SUNY Buffalo is acknowledged. The authors thank Dr. Nina Tyminska, Andrew Shamp, Zackary Falls, and Tyson Terpstra for their helpful comments.



REFERENCES

(1) Bartels, L. Tailoring Molecular Layers at Metal Surfaces. Nat. Chem. 2010, 2, 87−95. (2) Zhong, D.; Franke, J. H.; Podiyanachari, S. K.; Blömker, T.; Zhang, H.; Kehr, G.; Erker, G.; Fuchs, H.; Chi, L. Linear Alkane Polymerization on a Gold Surface. Science 2011, 334, 213−216. (3) Yang, B.; Björk, J.; Lin, H.; Zhang, X.; Zhang, H.; Li, Y.; Fan, J.; Li, Q.; Chi, L. Synthesis of Surface Covalent Organic Frameworks via Dimerization and Cyclotrimerization of Acetyls. J. Am. Chem. Soc. 2015, 137, 4904−4907. (4) Basagni, A.; Sedona, F.; Pignedoli, C. A.; Cattelan, M.; Nicolas, L.; Casarin, M.; Sambi, M. Molecules-Oligomers-Nanowires-Graphene Nanoribbons: A Bottum-Up Stepwise On-Surface Covalent Synthesis Preserving Long Range Order. J. Am. Chem. Soc. 2015, 137, 1802− 1808. (5) Gao, H. Y.; Held, P. A.; Knor, M.; Mück-Lichtenfeld, C.; Neugebauer, J.; Studer, A.; Fuchs, H. Decarboxylative Polymeriation of 2,6-Naphthalenedicarboxylic Acid at Surfaces. J. Am. Chem. Soc. 2014, 136, 9658−9663. (6) Sun, Q.; Zhang, C.; Cai, L.; Xie, L.; Tan, Q.; Xu, W. On-surface formation of two-dimensional polymer via direct C-H activation of metal phthalocyanine. Chem. Commun. 2015, 51, 2836−2839. (7) Shchyrba, A.; W ackerlin, C.; Nowakowski, J.; Nowakowska, S.; Bj ork, J.; Fatayer, S.; Girovsky, J.; Nijs, T.; Martens, S. C.; Kleibert, A.; 6639

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C

(29) Pal, A.; Vanka, K. Can silylenes rival transition metal systems in bond-strengthening [small pi]-back donation? A computational investigation. Chem. Commun. 2014, 50, 8522−8525. (30) Ruiz, V. G.; Liu, W.; Zojer, E.; Scheffler, M.; Tkatchenko, A. Density-Functional Theory with Screened van der Waals Interactions for the Modeling of Hybrid Inorganic-Organic Systems. Phys. Rev. Lett. 2012, 108, 146103. (31) Kunkel, D. A.; Hooper, J.; Simpson, S.; Rojas, G. A.; Ducharme, S.; Usher, T.; Zurek, E.; Enders, A. Proton Transfer in SurfaceStabilized Chiral and Polar Motifs of Croconic Acid. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 041402(R)(1−4). (32) Klimes, J.; Michaelides, A. Perspective: Advances and challenges in treating van der Waals dispersion forces in density functional theory. J. Chem. Phys. 2012, 137, 120901(1−12). (33) Miller, D.; Simpson, S.; Tyminska, N.; Zurek, E. Benzene Derivatives Adsorbed to the Ag(111) Surface: Binding Sites and Electronic Structure. J. Chem. Phys. 2015, 142, 101924. (34) Kunkel, D. A.; Hooper, J.; Simpson, S.; Miller, D. P.; Routaboul, L.; Braunstein, P.; Doudin, P.; Skomski, R.; Zurek, E.; Enders, A. Selfassembly of strongly dipolar molecules on metal surfaces. J. Chem. Phys. 2015, 142, 101921(1−9). (35) Fraxedas, J.; Garcia-Gil, S.; Monturet, S.; Lorente, N.; Fernandez-Torrente, I.; Franke, K. J.; Pascual, J. I.; Vollmer, A.; Blum, R. P.; Koch, N. Modulation of Surface Charge Transfer through Competing Long-Range Repulsive versus Short- Range Attractive Interactions. J. Phys. Chem. C 2011, 115, 18640−18648. (36) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; Fonseca Guerra, C.; van Gisbergen, S. J. A.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931−967. (37) Michaelides, A.; Ranea, V. A.; de Andres, P. L.; King, D. A. General Model for Water Monomer Adsorption on Close-Packed Transition and Noble Metal Surfaces. Phys. Rev. Lett. 2003, 90, 216102. (38) Mitoraj, M. P.; Michalak, A.; Ziegler, T. A Combined Charge and Energy Decomposition Scheme for Bond Analysis. J. Chem. Theory Comput. 2009, 5, 962−975. (39) Szatylowicz, H.; Krygowski, T. M.; Fonseca Guerra, C.; Bickelhaupt, F. M. Complexes of 4-substituted phenolates with HF and HCN: Energy Decomposition and Electronic Structure Analyses of Hydrogen Bonding. J. Comput. Chem. 2013, 34, 696−705. (40) Fonseca Guerra, C.; Zijlstra, H.; Paragi, G.; Bickelhaupt, F. M. Telomere Structure and Stability: Covalency in Hydrogen Bonds, Not Resonance Assisted, Causes Cooperativity in Guanine Quartets. Chem. - Eur. J. 2011, 17, 12612−12622. (41) Broclawik, E.; Zalucka, J.; Kozyra, P.; Mitoraj, M.; Datka, J. Formaldehyde activation by Cu(I) and Ag(I) sites in ZSM-5: ETSNOCV analysis of charge transfer processes. Catal. Today 2011, 169, 45−51. (42) Broclawik, E.; Zalucka, J.; Kozyra, P.; Mitoraj, M.; Datka, J. New Insights into Charge Flow Processes and Their Impact on the Activation of Ethene and Ethyne by Cu(I) and Ag(I) Sites in MFI. J. Phys. Chem. C 2010, 114, 9808−9816. (43) Rejmak, P.; Mitoraj, M.; Broclawik, E. Electronic view on ethene adsorption in Cu(I) exchanged zeolites. Phys. Chem. Chem. Phys. 2010, 12, 2321−2330. (44) Schleyer, P. R.; Maerker, C.; Dransfeld, A.; Jiao, H.; Hommes, N. J. R. Nucleus-Independent Chemical Shifts A Simple and Efficient Aromaticity Probe. J. Am. Chem. Soc. 1996, 118, 6317−6318. (45) Chen, Z.; Wannere, C. S.; Corminboeuf, C.; Puchta, R.; Schleyer, P. R. Nucleus-Independent Chemical Shifts (NICS) as an Aromaticity Criterion. Chem. Rev. 2005, 105, 3842−3888. (46) Fallah-Bagher-Shaidaei, H.; Wannere, C. S.; Corminboeuf, C.; Puchta, R.; Schleyer, P. v. R. Which NICS Aromaticity Index for Planar pi Rings is Best? Org. Lett. 2006, 8, 863−866. (47) Glassey, W. V.; Hoffmann, R. A Comparative Study of the CO/ M(111), M = Pt,Cu,Al Chemisorption Systems. J. Phys. Chem. B 2001, 105, 3245−3260. (48) Kresse, G.; Hafner, J. Ab initio molecular dynamics for liquid metals. Phys. Rev. B: Condens. Matter Mater. Phys. 1993, 47, 558−561.

et al. Controlling the Dimensionality of On-Surface Coordination Polymers via Endo- or Exoligation. J. Am. Chem. Soc. 2014, 136, 9355− 9363. (8) Lafferentz, L.; Eberhardt, V.; Dri, C.; Africh, C.; Comelli, G.; Esch, F.; Hecht, S.; Grill, L. Controlling on-surface polymerization by hierarchical and substrate-directed growth. Nat. Chem. 2012, 4, 215− 220. (9) Blyholder, G. Molecular Orbital View of Chemisorbed Carbon Monoxide. J. Phys. Chem. 1964, 68, 2772−2777. (10) Dewar, M. J. S. Bull. Soc. Chim. Fr. 1951, 18, C79. (11) Chatt, J.; Duncanson, L. A. Olefin co-ordination compounds. Part III. Infra-red spectra and structure: attempted preparation of acetylene complexes. J. Chem. Soc. 1953, 2939−2947. (12) Triguero, L.; Föhlisch, A.; Väterlein, P.; J, H.; Weinelt, M.; Pettersson, L. G. M.; Luo, Y.; Agren, H.; Nilsson, A. Direct Experimental Measurement of Donation/Back−Donation in Unsaturated Hydrogacron Bonding to Metals. J. Am. Chem. Soc. 2000, 122, 12310−12316. (13) Simpson, S.; Zurek, E. Substituted Benzene Derivatives on the Cu(111) Surface. J. Phys. Chem. C 2012, 116, 12636−12643. (14) Ö ström, H.; Nordlund, D.; Ogasawara, H.; Weiss, K.; Triguero, L.; Pettersson, L. G. M.; Nilsson, A. Geometric and chemical bonding of acetylene adsorbed on Cu(110). Surf. Sci. 2004, 565, 206−222. (15) Dyer, M. S.; Robin, A.; Haq, S.; Raval, R.; Persson, M.; Klimeš, J. Understanding the Interaction of the Porphyrin Macrocycle to Reactive Metal Substrates: Structure, Bonding, and Adatom Capture. ACS Nano 2011, 5, 1831−1838. (16) Rojas, G.; Simpson, S.; Chen, X.; Kunkel, D.; Nitz, J.; Xiao, J.; Dowben, P. A.; Zurek, E.; Enders, A. Surface State Engineering of Molecule-molecule Interactions. Phys. Chem. Chem. Phys. 2012, 14, 4971−4976. (17) Kunkel, D.; Simpson, S.; Nitz, J.; Rojas, G.; Zurek, E.; Routaboul, L.; Doudin, B.; Braunstein, P.; Dowben, P.; Eders, A. Dipole Driven Bonding Schemes of Quinoid Zwitterions on Surfaces. Chem. Commun. 2012, 48, 7143−7145. (18) Simpson, S.; Kunkel, D. A.; Hooper, J.; Nitz, J.; Dowben, P. A.; Routaboul, L.; Braunstein, P.; Doudin, B.; Enders, A.; Zurek, E. Coverage-Dependent Interactions at the Organics-Metal Interface: Quinonoid Zwitterions on Au(111). J. Phys. Chem. C 2013, 117, 16406−16415. (19) Cotton, F. A.; Wilkinson, G. Advanced Inorganic Chemistry; Wiley: 1988. (20) Hoffmann, R. A chemical and theoretical way to look at bonding on surfaces. Rev. Mod. Phys. 1988, 60, 601−628. (21) Crabtree, R. H. The Organometallic Chemistry of the Transition Metals; Wiley: 2005. (22) Churchill, M. R.; Mason, R. The Structural Chemistry of Organo-Transition Metal Complexes: Some Recent Developments. Adv. Organomet. Chem. 1967, 5, 93−135. (23) Valcarcel, A.; Clotet, A.; Ricart, J. M.; Delbecq, F.; Sautet, P. Comparative DFT study of the adsorption of 1,3-butadiene, 1-butene and 2-cis/trans-butenes on the Pt(111) and Pd(111) surfaces. Surf. Sci. 2004, 549, 121−133. (24) Yang, X. F.; Wang, A. Q.; Wang, Y. L.; Zhang, T.; Li, J. Unusual Selectivity of Gold Catalysts for Hydrogenation of 1,3-Butadiene toward cis-2-butene: A Joint Experimental and Theoretical Investigation. J. Phys. Chem. C 2010, 114, 3131−3139. (25) Braunschwig, H.; Damme, A.; Dewhurst, R. D.; Vargas, A. Bond-strengthening pi backdonation in a transition-metal pi-diborene complex. Nat. Chem. 2012, 5, 115−121. (26) Braunschweig, H.; Dewhurst, R. D. Single, Double, Triple Bonds and Chains: The Formation of Electron-Precise B-B Bonds. Angew. Chem., Int. Ed. 2013, 52, 3574−3583. (27) Himmel, H. Main group chemistry: Metal-reinforced bonding. Nat. Chem. 2013, 5, 88−89. (28) Hari Krishna Reddy, K.; Jemmis, E. D. Stabilization of diborane(4) by transition metal fragments and a novel metal to π Dewar-Chatt-Duncanson model of back donation. Dalton. Trans. 2013, 42, 10633−10639. 6640

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641

Article

The Journal of Physical Chemistry C (49) Blöchl, P. Projector augmented-wave method. Phys. Rev. B: Condens. Matter Mater. Phys. 1994, 50, 17953. (50) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (51) Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787−1799. (52) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. (53) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 38, 3098−3100. (54) Klimes, J.; Bowler, D. R.; Michaelides, A. Van der Waals density functionals applied to solids. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 195131(1−13). (55) Dion, M.; Rydberg, H.; Schröder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (56) Thonhauser, T.; Cooper, V. R.; Li, S.; Puzder, A.; Hyldgaard, P.; Langreth, D. C. Van der Waals density functional: Self-consistent potential and the nature of the van der Waals bond. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 76, 125112. (57) Román-Pérez, G.; Soler, J. M. Efficient Implementation of a van der Waals Density Functional: Application to Double-Wall Carbon Nanotubes. Phys. Rev. Lett. 2009, 103, 096102. (58) Lee, K.; Murray, E. D.; Kong, L.; Lundqvist, B. I.; Langreth, D. C. Higher-accuracy van der Waals density functional. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 081101(1−4). (59) Baerends, E. J.; Autschbach, J.; Bérces, A.; Bickelhaupt, F. M.; Bo, C.; Boerrigter, P. M.; Cavallo, L.; Chong, D. P.; Deng, L.; Dickson, R. M.; et al. ADF2010.01; http://www.scm.com. (60) Lenthe, E.; Baerends, E. J.; Snijders, J. G. Relativistic Regular Two-component Hamiltonians. J. Chem. Phys. 1993, 99, 4597−4610.

6641

DOI: 10.1021/acs.jpcc.6b00360 J. Phys. Chem. C 2016, 120, 6633−6641