Chemical insight from density functional modeling of molecular adsorption: Tracking the bonding and diffusion of anthracene derivatives on Cu(111) with molecular orbitals

Density functional theory (DFT) simulations have become a ubiquitous and indispensable tool in many areas of chemistry.¹ In surface science, they have been used to determine the adsorption geometry of organic molecules on metal, semiconductor, and insulator surfaces. For a vast array of adsorbate-substrate combinations, geometric optimization has revealed configurations that are consistent with experimental observations such as scanning tunneling microscopy (STM) imaging.² These simulations are based on minimization of the total energy of appropriately chosen super cells which contain the adsorbate as well as several layers of substrate lateral unit cells. This approach has led to the understanding of surface pattern formation in a plethora of adsorption systems;^3,4 likewise, it has permitted faithful simulation of molecular motion at surfaces.^4,5 The overall success of this method is beyond doubt. Yet, from a chemical perspective, this approach is somewhat unsatisfactory, as interpretation of its results, i.e., charge-density distributions, and their reconciliation with a molecular-orbital (MO) approach to chemical understanding is difficult and cumbersome and, hence, all too often not even attempted. Here, we present a facile and very informative way of analyzing DFT data on molecular adsorption at metal surfaces that provide a direct bridge to a more chemical and MO-oriented understanding of the underlying interactions.⁶ Appendix  C discusses the code implementation of our approach so that it can be applied to general systems.

While our approach does not generate novel or better data per se, we submit that it improves chemical intuition about adsorption as well as surface diffusion processes and, thus, may lead to deeper insights into surface chemistry. In particular, we demonstrate a way to obtain a MO picture—including approximations to MO diagrams—by utilizing the projections of the Kohn-Sham (KS) orbitals onto the partial waves of each atom, which is generated by typical DFT codes using projector augmented waves (PAW).^7,8 This permits highlighting the impact of individual functional groups in the adsorbate-substrate interaction, directly relating DFT-analysis to conventional chemical intuition; in fact, application of KS orbitals to a MO analysis is perfectly reasonable,^9,10 if not more appropriate than some of the more traditional orbitals (e.g., Hartree-Fock) that are typically applied to a qualitative MO theory picture. Indeed, our approach may give pause to those who view KS orbitals as purely formal entities devoid of intrinsic physical meaning, while not surprising the many practitioners who recognize how good DFT-based orbital descriptions can be.

As examples, we apply our method to a family of anthracene derivatives adsorbed on the Cu(111) surface: anthraquinone (AQ), dithioanthracene (DTA), and diselenoanthracene (DSeA), as well as anthracene (A) itself (a fairly rigid linear aggregate of three fused benzene rings, Fig. 1(a)). AQ, DTA, and DSeA are derived from anthracene by replacing the hydrogen atoms at the central 9,10-positions with the elements O, S, and Se, respectively, which are the upper three members of the chalcogen group of the periodic table. Because of the chalcogens’ comparatively strong interaction with the substrate, we refer to them as linkers in this work. While our DTA and DSeA precursors have extra functional groups attached to their linkers (in order to prevent polymerization), for the experiments discussed here we have removed those groups by annealing after deposition on Cu(111). Thus, for the experimental studies of these systems to which we have applied our MO analysis,^11,12 all three species have the same basic structure.

This class of organic molecules deposited onto Cu(111) is of particular interest because its members have shown some promise as prototypical molecular machinery: they are able to diffuse uniaxially along the [110] direction of the Cu surface and, in the case of AQ, even act as a cargo-carrying agent.^11,12 Both AQ and DTA are also capable of creating long-range ordered structures on the Cu(111) surface (hydrogen-bonded chains¹³ in the cases of both DTA and AQ, and a giant honeycomb network^14,15 in the case of AQ). Anthracene derivatives, in general, have also found applications in photovoltaic devices¹⁶ and electronics.¹⁷ 

In this work, we first examine the geometric and energetic details of adsorption for these three species as well as unsubstituted anthracene. We then reexamine adsorption from the point of view of the molecules’ KS orbitals, which we show can readily be partitioned into contributions from KS states of the anthracene backbone and contributions from KS states derived from the linkers. Thus, we can identify which of these states participate in hybridization upon adsorption (i.e., which molecule states are most involved in bonding to the substrate). As a consequence, we can systematically (in a chemical sense) analyze differences in the adsorption-related orbitals and bond configurations. Variations within this set offer a way to understand the origins of each molecule’s diffusion barrier, as extracted from experimental STM data for diffusion of each species; some aspects of anthracene, AQ, and DTA (but not DSeA) have already been published,^11,12 and others are readily found in this study via DFT including (experimental/computational) diffusion barriers of 30 meV/72 meV, 130 meV/135 meV, [too low to be resolved; cf. Sec. VII]/79 meV for AQ, DTA, and DSeA, respectively.

The principal finding of this study is that differences in diffusion barriers of the species can be correlated to different linker-derived states hybridizing to the underlying Cu(111) substrate. Side-by-side comparison of MO diagrams of each adsorbed system permits analysis of how each orbital contributes to the overall diffusion barrier. In particular, the chalcogen linkers of DSeA are found to bind in a fundamentally different way from those of AQ and DTA, suggesting an explanation for the low DSeA diffusion barrier observed by STM (see below).

Our DFT calculations use the VASP code¹⁸ (version 5.2.12) incorporating non-local van der Waals (vdW) interactions in the vdW-DF1 approximation¹⁹ with the optB86b-vdW^20,21 exchange–correlation functional and a PAW basis.^7,22 (The Perdew-Burke-Ernzerhof Generalized Gradient Approximation (PBE-GGA)²³ functional was also used to measure binding energies for comparison.) The optB86b-vdW implementation was highly rated at the time when these calculations were begun²¹ and has continued to test well for height and binding energies in adsorption of large organic molecules on Cu(111).²⁴ In the meantime, the family of vdW-DF functionals has grown, and presumably one could utilize its other members such as the vdW-DF2 formulation,²⁵ which is excellently suited to this problem,²⁴ especially once Hamada’s very recently revised form of the Becke 86 exchange²⁶ is included. Another option is the non-empirical vdW-DF-cx formulation,^27,28 which is very promising and is currently being benchmarked. We refer readers to an extensive review about to appear.²⁹ 

Calculations were performed for each molecule deposited onto an 8-layer Cu slab, with each layer having 6 × 5 atoms lateral periodicity. Dipole corrections were employed to cancel the effects of any induced dipoles: the charge distribution center was chosen to be the geometric center of the atomic coordinates, energetic corrections were calculated for a dipole along the z-direction of Fig. 1(b) (perpendicular to the slab), and standard corrections³⁰ to the periodic potential that cancels out the net dipole between super cells were employed.³¹ All results were optimized so that the remaining forces were less than 0.03 eV/Å Because of the large size of the unit cell and the ensuing small volume of the first Brillouin zone, we sampled k-space exclusively at the $Γ ¯$ point. Plane-wave and augmentation-charge energy cutoffs of 400 and 700 eV, respectively, were used. A detailed explanation of our MO analysis is presented in Sec. IV.

Before developing a systematic molecular-orbital picture of the variations of the substrate-linker-backbone or of the adsorbate-substrate interactions in each adsorption system, we first describe the overall adsorption geometry and associated binding energies. Figures 1(a) and 1(b) show top and side views, respectively, of what we consider to be the minimum-energy adsorption configuration assumed by anthracene derivatives on the Cu(111) surface. One linker lies at an atop site while the other lies at or near a fcc hollow site; the rings of the anthracene backbone are centered above hcp hollows. This determination was made by noting that the long axis of the anthracene backbone points along the close-packed [110] direction in STM images^11,12 and then ascertaining the corresponding fully relaxed configuration by DFT.

The binding-energy graph of Fig. 1(c) shows that each anthracene-derived molecule is stable on the surface with a binding energy that increases as the linker species proceeds down the periodic table from H via O to Se. The importance of incorporating a vdW-DF functional has become well recognized.^28,32 Calculations without the vdW-DF corrections (PBE-GGA) yield a low binding energy for AQ—on the order of a few hundred meV—which is inconsistent with experiment in which AQ sticks on Cu(111) up to and above room temperature.¹⁴ This result is in line with other non-vdW DFT studies of aromatic systems on metal surfaces in which calculated binding energies are found to be underestimated when compared to experimentally determined ones and have been linked to a lack of proper accounting of the van der Waals forces,^4,33 as well as more recent studies³⁴ in which inclusion of van der Waals interactions has shown good agreement with experiment. The difference in calculated binding energy between the optB86b-vdW functional and the PBE-GGA functional is 1.5 eV for A, while for each of its derivatives the difference is approximately 2 eV (1.99, 2.11, and 2.08 eV for AQ, DTA, and DSeA, respectively; cf. Fig. 1); these differences can be viewed as indicative of the total dispersive interaction for each of the species considered, suggesting that for each of the three derivatives, the dispersive interaction is similar and clearly a necessary ingredient in any description of their behavior.

Consideration of the adsorption geometries of the molecule and substrate sheds further light on the trend in binding energies. Since PBE-GGA has difficulty in binding AQ, we present only the results for the optB86b-vdW functional here; we also find essentially similar results using the optB88-vdW functional. The linker atoms of each anthracene derivative are pulled close to the substrate (compare the linker z-coordinates L_1z and L_2z to the central ring heights C_z of each molecule in Fig. 1(c)), while the anthracene backbone bends upward in a V-shape, as can be seen in the side view of Fig. 1(b), in which the z direction has been exaggerated for clarity. The values of C_z in Fig. 1(c) show that unmodified anthracene, in contrast, comes to rest more than 3 Å above the surface. Thus, the V-shape of the derivative molecules, which appears as a set of 2 lobes in STM images,^11,14 occurs because the linkers pull the central ring closer to the substrate, with the ends of the anthracene backbone then bending upward. The relative degree of bending for each molecule is indicated by the quantity T_z − C_z, the height of the terminal carbon atoms above the central ring; for the unmodified anthracene, T_z − C_z vanishes (i.e., the anthracene is flat), in agreement with experiment.³⁵ 

In the presence of a linker to the substrate, the central ring height tends to increase with increasing linker atomic radius (i.e., from O to Se). The linker size may also help to explain the same trend in binding energy: larger linkers require less forcing of the anthracene backbone to the substrate, contributing to higher net adsorption energies and a smaller height difference between T_z and C_z, corresponding to less bending of the molecule (bottom right axis of Fig. 1(c)). The substrate distortion, S, measured here as the z-coordinate of the most uplifted Cu atom, also appears to mirror this trend of decreasing distortion going from AQ to DSeA; however, DTA shows more substrate distortion than AQ; a possible explanation for the increased distortion is presented in Sec. VII.

In order to elucidate the differing impacts of the chalcogen linkers, we need to identify which electronic states dominate the substrate interactions and how these states shift in energy and real space as the linkers are varied. To this end, we take the KS states of each system and create the site-projected overlap integrals between pairs of states (one from each system): if the squared magnitudes of their projections onto one another are above a threshold within the PAW spheres of the atoms unchanged by substitution, then they are considered to be either the same state or, at the very least, derived from one another. Here, PAW spheres, defined automatically in typical PAW codes, refer to the regions surrounding each atom over which the projectors are defined (those regions where electronic wave-functions would typically oscillate rapidly).

For example, if we wish to relate unsubstituted A to AQ in vacuum, we perform DFT calculations on these two related systems: AQ isolated in vacuum and A isolated in vacuum. Then, since the hydrogen and carbon atoms of A remain unaltered in the substitution that yields AQ (with the exception of the 9,10-H atoms which are replaced by O atoms), the projection of the state $| ψ n AQ 〉$ from AQ onto the state $〈 ψ m A |$ is given by

where the first sum is over the atoms, a, that remain unchanged (in this example, the carbon and most of the hydrogen atoms), the second sum is over the partial wave indices, i, centered at each of these atoms, and $| p ̃ i a 〉$ are the PAW projector functions defined by Blöchl.⁷ Note that the term on the left hand side of Eq. (1) is not the full projection of $| ψ n AQ 〉$ onto $〈 ψ m A |$ but rather the projection within the PAW spheres; i.e., the sum over the density matrix of these PAW functions is not unity. We, therefore, replace the vertical line by a vertical ellipsis; this allows readers to benefit from the intuitive convenience of imagining this quantity as a projection while providing a constant reminder that it is not a full projection. Further discussion of the approximate nature of Eq. (1) is given in Appendix  A. Other aspects of the notation presented here are the same as adopted in GPAW^36,37 and by Blöchl⁸ except that the index corresponding to the atomic site is written explicitly as a superscript a, rather than being incorporated into the subscript i. Because the quantities $〈 ψ ̃ m A | p ̃ i a 〉$ and $〈 p ̃ i a | ψ ̃ n AQ 〉$ are already calculated as part of any PAW code performing DFT on the systems of anthracene and AQ, respectively, they are readily available as output from such calculations.³⁸ 

Since the projector functions are zero outside of the PAW spheres, they do not form a complete set over all of space, necessitating that the resulting projections be normalized in a way that correctly indicates the degree to which the two states being compared can be considered as well represented by one another (i.e., they overlap sufficiently). We find reasonable results if we construct an ad hoc normalization constant, d_mn, equal to the maximum of $|〈 ψ m A ⋮ ψ j AQ 〉 | 2$ and $|〈 ψ j A ⋮ ψ n AQ 〉 | 2$⁠, where j runs over all of the KS states of AQ and anthracene, respectively, and use it to define a matrix representation of the MO diagram

so that if $ψ m A$ and $ψ n A Q$ are identical, MO_mn is unity.

For those pairs of states that are similar in spatial extent in the systems of AQ and A (one state from each system), only the near diagonal elements of the MO matrix are non-zero (since the matrix is not square, the diagonal is not precisely defined), whereas states affected by the insertion of the functional groups are represented by rows with off-diagonal elements. The hybridization partners can then be read from the position of the non-zero elements in the row. The MO matrix can consequently be compactly represented as a MO diagram. To simplify the diagrams, we additionally set values of MO_mn to 0 if they fall below some threshold. The threshold value determines how much detail to be displayed in the MO diagram and, thus, depends on the system studied and the objective at hand. In this work, we chose values of 0.1 and 0.35 (out of a maximum of 1) for systems in vacuum and for adsorption systems, respectively, so as to generate diagrams that are visually clear yet omitting none of the important features of the electronic structure. Examples of these diagrams for AQ, DTA, and DSeA are in Figs. 2 and 4, the details of which are explained in Secs. V and VII, respectively. Each column of line segments shown represents a different MO matrix linking the states of the system on the left to those on the right.

Figure 2 shows the relationships between KS orbitals of the A precursor and each of its derivatives, which in Sec. VII are found to hybridize extensively with the Cu(111) surface upon adsorption (i.e., the MO matrix has many off-diagonal elements). Fig. 2(a) was generated using the implementation described in Sec. IV, and as such serves as a first test case of relatively simple systems whose orbital relationships can easily be verified by visual inspection (see the surfaces of constant KS orbital charge density of Fig. 2(b)). The first column of Fig. 2(a) shows the density of states (DOS) of A. The second column shows the DOS of AQ. Between these two columns, lines connect the eigenvalues of KS states of A to those states of AQ that derive from them. Similarly, the DOS of each anthracene derivative is shown in its respective column, and the related states in each can be traced all the way (from right to left) back to their A precursor states. Those states that we consider to have some importance are color-coded according to Fig. 2(b). The order of color assignment for each state of each system is based on which states between systems are comparable (i.e., which states are linked together). For example, in Fig. 2(a), since the LUMO of A clearly connects to the LUMO + 1 of AQ, both states are colored purple; geometric comparison between the LUMO of A and the LUMO + 1 of AQ in Fig. 2(b) shows that this assignment is correct. We use red, orange, and yellow for new states contributed from the linkers that do not exist in the A precursor (whose states are represented by bluish colors). While the HOMO-4 of AQ (red, second column) derives partially from the A HOMO-4, it is color-coded as a linker-derived state because the equivalent state in DTA and DSeA shows little relation to any state from A, as can be seen in the HOMO-1 states of Fig. 2(b) outlined in red. We treat the HOMO in the derived systems as being degenerate and so assign it two colors (yellow and orange). The repetition of the color blue for the LUMO and HOMO-1 of AQ (HOMO-2 for DTA and DSeA) indicates that both of these states are derived from the HOMO of A.

As already noted, Fig. 2(a) could have been created by direct inspection of the shapes of the KS orbital densities of Fig. 2(b). This exercise is intended as a proof of concept that the procedure described in Sec. IV is capable of making the correct assignments, as well as establishing the similarities between all three anthracene derivatives and their relationship to the A precursor. Figure 2(b) also makes it clear that the A HOMO-derived states (boxed in blue) are hybridizations with the p_z orbitals of the linkers, creating a π bonding and anti-bonding pair, with the implication that these states, while assigned to the A precursor, should also have properties related to the linkers. For the class of systems addressed experimentally in this work, this analysis establishes that our PAW-based method can correctly obtain the overlaps despite its inherent approximations. Application of this method (as well as its further applications described in Secs. VII and VIII) to systems containing other groups of atoms will await further experimental work to provide a direct validation.

Having introduced the molecular orbitals of anthracene and its derivatives in vacuum, we now motivate application of this same procedure to these species adsorbed on Cu(111) by first performing a real-space analysis of charge transfer. This approach illuminates fundamental aspects of the adsorption process and at the same time sets the stage for the study of the diffusive behavior of the anthracene derivatives.

Figure 3(a) introduces the three energy minima of the potential energy surface (PES) for diffusion of each of the anthracene derivatives along its long axis, the preferred diffusion direction at least for AQ and DTA. As described at length in previous work,^11,12,39 the admolecules diffuse by initially stepping one linker into a hollow site via a shallow barrier, followed by the main barrier when the other linker swings around to place the molecule at a fcc site in the transition state. This is followed by a final movement of the first linker from its hollow site across the same shallow barrier to return to a hcp site equivalent to the initial one. We note that the tilted intermediate minimum adopted by AQ differs gradually from those adopted by DTA and DSeA. The hcp and fcc configurations are named so because the centers of the aromatic rings fall above hcp and fcc hollows, respectively, while the tilt configuration lies between the two, centered over a bridge site. The fcc configuration is calculated to be the highest in energy along the path of minimum energy, and as such can be considered to be the transition state (at the top of the diffusion barrier), while the hcp configuration is the ground-state configuration, as noted earlier.

Figure 3(b) shows differential charge density iso-surfaces for each system in each configuration, calculated by subtracting the charge densities of the Cu slab and bent molecule (as independent systems calculated separately in vacuum) from the charge density of the adsorbed (combined) system. Regions of negative charge density difference then correspond to regions from which electron density was lost, while those of positive charge density difference correspond to regions that gained electron density. The colors applied to each image (see bottom inset of Fig. 3(b)) correspond to a Hirshfeld partitioning⁴⁰ of the charge density into parts corresponding to the aromatic rings (the backbone), the linkers, and the Cu substrate. From Fig. 3(b), it is apparent that in moving from the hcp to the fcc configuration (following one row from left to right), the regions defined by the iso-surfaces do not change much in size or shape beyond deforming around the motions of the linkers and the anthracene backbone, suggesting that in order to understand the differences in diffusion of these species, it is sufficient to characterize the differences in how each molecule bonds to the substrate in just one of the three configurations.

Qualitatively, the geometries of the iso-surfaces (Fig. 3(b)) resemble the vacuum KS orbitals of Fig. 2(b). This correspondence is most striking when illustrated as in Fig. 3(c), which shows a side by side comparison of the charge density gain of DSeA upon adsorption (left) to the LUMO of DSeA in vacuum (right). The similarities between the two strongly suggest that the LUMO of DSeA is being filled upon adsorption to the Cu(111) surface. Similar behavior is also observed for DTA and AQ.

We now apply the MO analysis of Sec. IV to the adsorbed derivative systems; these systems exhibit a continuum of states in principle and on the order of thousands of KS states in practice (i.e., in VASP calculations), making the automated analysis described in Sec. IV indispensable. Figure 4 shows comparisons between the adsorption (in the hcp adsite configuration) of (a) AQ and DTA and of (b) AQ and DSeA. It is arranged in a manner similar to Fig. 2(a), with the central DOS graphs in each panel being the projected density of states (PDOS) onto the adsorbate atoms (i.e., excluding PDOS of the underlying Cu) of the adsorbed systems. The orbital color-codes introduced in Fig. 2 are retained and repeated in Fig. 4(c). The first two columns of Fig. 4(a) are repeated from Fig. 2(a) (showing how the states of AQ are derived from those of A). The third column illustrates how the KS orbitals of an isolated molecule of AQ hybridize with the underlying Cu to create the states of the adsorbed AQ system. The three right columns repeat the same analysis for DTA, but in reverse (right to left) order. This setup allows for a direct comparison of the adsorbed-system KS orbitals of AQ and DTA in the center column. Fig. 4(b) repeats the same structure, but now for direct comparison of adsorbed AQ (left) to adsorbed DSeA (right).

Tracing the path of the blue LUMO state in the third and fourth columns of both Figs. 4(a) and 4(b) shows that the LUMO of each derivative system is indeed filled upon adsorption (having moved from above the Fermi level in the vacuum system to below the Fermi level in the adsorbed system), as was independently verified by inspection in Fig. 3(c) of Sec. VI. From the differential charge density iso-surfaces of Fig. 3, it was not so obvious that the LUMO + 1 (purple) is also partially filled in each system. Additionally, this MO diagram analysis makes it easy to scrutinize the filled states. Their contributions to adsorption cannot be so easily gleaned from a differential charge density plot because they are not gaining or losing electrons to the extent of the LUMO (which goes from empty to nearly full as determined by calculating the overlap of the vacuum system orbital with its hybridized occupied orbitals in the adsorbed system). The filling of the LUMO and partial filling of the LUMO + 1 suggest a rather substantial electron transfer from the substrate to the molecule. However, this is offset by electron transfer from energetically lower-lying molecular states (as evidenced by the electron loss indicated with green and blue in Fig. 3(b)) in a manner analogous to the π bonding and back-bonding of CO to transition metals as in the Dewar–Chatt–Duncanson (DCD)/Blyholder model;⁴¹ such to-and-fro charge transfer has also been observed for various organic molecules on transition metal surfaces.⁴² 

Scrutiny of the third columns from left and right of Fig. 4(a) shows that the bonding of AQ and of DTA to the surface is quite comparable, particularly for occupied states. While the precise KS eigenvalues vary slightly between the two, they have a near one-to-one correspondence; what does appear to be an important difference, however, is the tendency of the HOMO and HOMO-1 (orange and red states, respectively) of DTA to hybridize to comparatively lower KS energies as highlighted by the dashed ovals of Fig. 4(a). The increased substrate distortion associated with DTA when compared to AQ (see Sec. III) is likely associated with this stronger binding of DTA’s linkers.

Fig. 4(b) shows that DSeA hybridizes quite differently from the other two derivatives, as is evident when comparing adsorbed AQ to adsorbed DSeA (and by extension, when comparing DTA to DSeA). In particular, DSeA shows significantly more hybridization of its LUMO, HOMO-1, HOMO-2, and HOMO-4 (blue, red, blue again, and cyan, respectively) in the higher range of −2 eV to 0 eV (circled). Furthermore, DSeA shows no signs of bonding of its HOMO to lower energy states (below −3 eV), whereas AQ and DTA do; these differences are highlighted by the dashed ovals of Fig. 4(b).

Here, we report that we performed experiments with DSeA on Cu(111). Figure 5 shows an image of linear molecular aggregates quite similar to those observed for DTA.¹³ Most striking, however, is a 2D gas of DSeA molecules (e.g., blurry/noisy regions in Fig. 5(a)) surrounding them that could not be frozen out at the lowest temperature of our instrument (15 K), in stark contrast to DTA and AQ whose motions are frozen out even at higher temperatures.^11,12 We observed a number of aggregation structures, which will be reported elsewhere. This finding suggests a significantly lower diffusion barrier than those previously reported for AQ¹² and DTA¹¹ and related molecules.⁴³ 

Relying on the KS energies rather than on more physical energies, the preceding description offers what is essentially a qualitative assessment. Consequently, in Sec. VIII, we present a quantitative approach to interpreting how these hybridizations affect diffusion from which it will become apparent that the differences highlighted by green circles in Figs. 4(a) and 4(b) are indicative of DTA’s comparatively higher diffusion barrier and DSeA’s comparatively lower barrier.

The links illustrated in Fig. 4 provide a way to trace back properties of the adsorbed system to the MOs of the vacuum system from which they derive, both visually as presented in Sec. VII and, as will be demonstrated, in terms of energetic contributions to the diffusion barrier. Here, we investigate which MOs attract a given anthracene derivative to the energetic minimum adsite (the hcp configuration) versus which MOs attract the molecule to the transition state adsite (the fcc configuration): it is the competition between these two bonds that determines the diffusion landscape. The KS energies of the adsorbed system represent non-interacting single-particle energies and, therefore, an inappropriate way to determine the contributions that their associated orbitals make to the diffusion barrier because they do not completely account for the energy of the fully interacting system (i.e., the sum of the KS energies does not yield the ground state energy of the system). It is, however, possible to assign energies to the MOs that do represent physical energies by making use of the Hellmann-Feynman theorem.⁴⁴ This is accomplished by assigning partitions of the total charge density of the adsorbed system to each vacuum-system MO according to the links between the MOs of Fig. 4. For a given MO, the energy it contributes when the atoms of the system move a small distance out of their preferred configurations and towards the transition state is simply the dot-product of the Hellmann-Feynman force exerted by its associated charge density partition on each nucleus and the displacement vector traveled. Appendix  B describes the formal underpinning of this analysis in detail. By integrating over the entire path of motion from the energetic minimum configuration to the transition state (from the hcp adsite to the fcc adsite in this case), the total energetic contribution of each molecular orbital is calculated, as described in Appendix  C; the results are presented in Fig. 6.

The magnitudes of the bars shown in Fig. 6 represent the relative energetic attraction that each molecule feels either towards the hcp adsite (bars pointing left) or towards the fcc adsite (bars pointing right) as contributed by each of its MO-derived charge distributions. According to this analysis, DSeA only feels attraction to the hcp adsite through its LUMO + 1 and HOMO-3 orbitals, while all others are acting to push it toward the fcc adsite. In comparison, AQ and DTA have half of their orbitals attracting them to each adsite, with DTA showing a strong attraction from both its HOMO and HOMO-1 (green oval in Fig. 6(b)) states toward the hcp adsite (the energetic minimum), resulting in the pronounced diffusion barrier observed experimentally; this observation agrees with the green ovals of Fig. 4(a), which highlight how at the hcp site the hybridized HOMO and HOMO-1 of DTA have lower KS energies than those of AQ and DSeA.

The tendency for DSeA’s MOs to favor the fcc configuration seems consistent with the observation that the HOMO-1, HOMO-2, and HOMO-4 of DSeA hybridize more in the higher KS energy range than those of AQ and DTA, as highlighted in the green ovals of Fig. 4(b). The same orbitals have been circled in Fig. 6(c), illustrating that they do indeed result in DSeA favoring the hcp configuration less than AQ and DTA do; if these orbitals behaved the same as they do in AQ and DTA, then DSeA would be more strongly attracted to the hcp site, and its overall diffusion barrier would be larger, more like the other two anthracene derivatives.

We also note the overall chemical similarity of the three linker atoms, which we attribute to their common chalcogen characteristics: their LUMOs have similar extent and fill in a closely related fashion upon adsorption; their preference for initial vs. transition state shows congruence (their bars point to the same side in the dashed gray box of Fig. 6). In combination, these observations have the important implication that to understand differences between these otherwise similar species, it is necessary to look beyond the frontier orbitals and, at the very least, to consider the behavior of other orbitals that are strongly involved in hybridization with the underlying substrate (i.e., LUMO + 1 through HOMO-4 in this study).

In conclusion, we have demonstrated a way to extract MO diagrams from typical DFT codes employing PAW bases for a variety of related systems that enables facile side-by-side comparison of their chemical properties within a framework based on MO theory. Using this technique, we were able to trace the orbitals of anthracene derivatives from their beginnings as anthracene orbitals and linker orbitals to their final configurations upon adsorption of each derivative. Based on the knowledge of the origins of each of the adsorption systems’ orbitals, we could partition the diffusion barrier energy into the separate/individual contributions from each vacuum-system MO. This theoretical analysis is complemented by experimental findings of rapid DSeA diffusion on Cu(111) even at very low temperatures, i.e., temperatures at which the motion of DTA and AQ were already frozen out.

We gratefully acknowledge support by DOE DE-FG02-07ER15842, as well as joint NSF Grant Nos. CHE 13-06969 (L.B.) and CHE 13-05892 (T.L.E.). This research used resources of the National Energy Research Scientific Computing Center, which is supported by DOE DE-AC02 05CH11231.

In this work, as described in Sec. IV, we utilize the quantities $〈 p ̃ i a | ψ ̃ n 〉$⁠, which are the auxiliary function projected onto a projector function in the language of Blöchl^7,8 or the projector overlaps in the language of GPAW.^36,37 They are readily available as output from the popular DFT code VASP.¹⁸ However, they constitute an approximation because they are defined solely within the PAW spheres centered about each atom. This raises concerns relating to the treatment of charge density outside these spheres and the overlap of PAW spheres. The former concern, though, is mitigated by the fact that the PAW procedure updates the description of core levels in response to the hybridization that occurs in the chemically active regions.

To address the second concern, we checked how the radii of these cutoff spheres compare to the typical inter-atomic distances in the systems studied. The PAW sphere radii (as defined by the RCORE values in VASP) for H, C, O, S, Se, and Cu, respectively, are (in Å) 0.58, 0.79, 0.80, 1.0, 1.1, and 1.2. These radii are sufficiently large that the sum of nearest neighbors’ values can be greater than their separation (by of order 0.1 Å). Kresse and Joubert²² and Furthmüller et al.⁴⁵ note that VASP calculations yield accurate results due to their particular construction of the projector functions despite these overlaps. We also find that plotting of the pseudo-wavefunctions from the POTCAR file reveals cutoffs not reported explicitly in the output that is, in fact, small enough to avoid overlap between nearest neighbors, suggesting that the projector functions implemented in VASP may indeed not overlap. The results in Sec. V further corroborate the earlier finding⁵ that the PAW projections allow, at the very least, qualitatively correct identification of which states of a substituted system are derived from a precursor system.

We note that the GPAW³⁷ code offers two techniques for calculating overlap between KS orbitals of two related systems: (1) a “projector-pseudo-wavefunction overlap” in which the projection of a full wavefunction onto an atomic orbital $ϕ i a$ is approximately given by $〈 ϕ i a | ψ n 〉≈〈 p ̃ i a | ψ ̃ n 〉$⁴⁶ and (2) a less efficient but more accurate grid-based “all electron (AE) wavefunction overlap,” which is available in GPAW, but not necessarily in other codes. As already noted, we use an approach related to the first method.

The KS orbitals are considered to be an auxiliary fictitious formal system used as a tool to solve the many-electron ground state problem in DFT. However, it is already common practice to treat KS orbitals as though they represent some real physical quantity (e.g., they are routinely used to generate DOS or PDOS plots). The main problem with these orbitals is that, because they represent the contrived non-interacting single-particle electrons, their energy eigenvalues are not, strictly speaking, the energy eigenvalues of actual particles (or quasi-particles). However, Stowasser and Hoffmann⁹ noted that the energy eigenvalues can often be rescaled to fit, for example, the Hartree-Fock orbital eigenvalues. Furthermore, a deeper interpretation of their physical meaning lies in the charge density which they can be considered to partition.

In MO theory, particular orbitals are typically characterized as “bonding,” “anti-bonding,” or “non-bonding,” the determination of which is often made on a geometric basis (viz., whether the charge distribution of a particular orbital helps to hold 2 nuclei together, pushes them apart, or does nothing significant, respectively). For complex systems such as those studied in our work, this definition presents a problem because there are many nuclei. The designation of bonding or anti-bonding more generally indicates whether the orbital promotes or hinders, respectively, some specified change to the system, in particular, here, adsorption on a specific kind of site. Thus, one orbital that binds a molecule to a substrate might prefer to bind the molecule to one adsite is said to be bonding with respect to that site but could be antibonding with respect to another; likewise, a different orbital might be bonding with respect to the second site, etc. (In other words, bonding and antibonding are synonymous with promoting or hindering, respectively, for some system modification.) We seek to analyze the PES associated with motion of each of the nuclei (typically easily calculated with DFT, though famously problematic for some systems⁴⁷) to determine the contribution of each KS orbital. Thus, a diffusion barrier (see Fig. 7) can be decomposed into its energetic contributions from each KS orbital with orbitals that increase that barrier considered as bonding with respect to the preferred adsite (A in Fig. 7) and with those orbitals diminishing the barrier considered as anti-bonding with respect to the preferred adsite (or, equivalently, bonding with respect to the transition state (B)). Hence, the quantitative energetic role of each KS state (or, more generally, any chosen partitioning of the charge density) is manifested. We show below how to do this decomposition and classification.

Based on the Hellmann-Feynman theorem,⁴⁴ the ground-state charge density can be evaluated directly as an electrostatic charge distribution in order to calculate the force on each nucleus (or ion in a frozen-core approximation) in a given system. In this way, one can, in principle, calculate the energy difference between, for example, two different adsorption configurations (call them A and B) of a molecule undergoing diffusion (i.e., we are considering the PES at two points A and B). One would do this by dividing up the path from A to B into differential changes of the positions of the nuclei and then calculating the work done in forcing the ions to move from configuration A to B

where the index i runs over the combination of each of the nuclear coordinates as well as Cartesian components, x_i (e.g., one might index such that x₀, x₁, x₂ are the x,y,z coordinates, respectively, of the 1st nucleus, x₃, x₄, x₅ are the x,y,z coordinates of the 2nd nucleus, etc.). The force component on the nucleus at each such point is F_xi, which is then split into a contribution from the ground-state electron density $F x i elec n 0$ evaluated at each differential configuration and a contribution from the repulsion due to the surrounding nuclei $F x i nuc$⁠. The energy difference between A and B could, of course, be determined by just taking the difference in the ground-state energy calculated by DFT for the system in the two configurations. However, the first approach has the advantage that the energy has been written in a way that can be decomposed by the KS orbitals; since $F x i elec$ is linear with respect to partitioning of the charge density, we can define the force on each ion due to the occupied KS orbital |ε〉 as $F x i ε ≡ F x i elec n ε$⁠, where $n ε ( r ⇀ ) =|〈 r ⇀ ∣ε〉 | 2$⁠; explicitly,

Alternatively and advantageously, we can consider the force associated with an orbital |e〉 of the related antecedent molecule-in-vacuum system: $F x i e , ε ≡〈ε|e〉〈e|ε〉 F x i ε$ so that $F x i ε = ∑ e F x i e , ε$⁠. The advantage is that |e〉 does not change as the adsorption system goes from A to B, so that substitution back into Eq. (B1) yields quantities (denoted $Δ E A → B e$⁠) very similar to what we desire

where the expression for $Δ E A → B e$ is given by the term in braces above. This is already a satisfactory and physically interpretable partitioning of the barrier energy (if A is taken to be the optimum configuration and B the transition configuration) into parts due to the |e〉 states and one large contribution due to the nuclear repulsions. However, to be more faithful to the original concept of bonding/anti-bonding, our definition of the effective energy that a KS orbital contributes to the barrier includes the work done by each orbital in opposing the repulsions of the nuclei (given N occupied states):

and

where $Δ E ̃ A → B e$ is again given by the term in braces in Eq. (B3), and the use of the tilde indicates effective force components and energies. Now for each antecedent molecule state |e〉 (LUMO, HOMO, etc.), an energy contribution to the barrier can be defined to create a picture similar to Fig. 7 (see Fig. 6), yielding at the same time a quantitative interpretation of the MO type analysis and the KS orbitals themselves: the calculated quantity $Δ E ̃ A → B e$ allows us to tabulate the change in energy due to the electrostatic forces imparted by the subset of KS states stemming from precursor state |e〉. For example, the upward pointing green arrow of Fig. 7 labeled as ΔE LUMO indicates a positive value of $Δ E ̃ A → B LUMO$⁠, implying that when the LUMO hybridizes with the underlying substrate, the resulting states bind it more strongly to site A than to site B.

In this section, we provide a brief description of how the technique of Appendix  B is implemented without modifying the source code of VASP (i.e., an implementation that only requires processing the output files of VASP calculations). Most of the necessary processings are performed by code that we have made available.⁴⁸ This code includes a Mathematica® (Ref. 49) notebook (GenMODiagrams.nb) that is used to (1) generate the appropriate VASP input files, (2) generate shell scripts (unix/linux bash scripts) for executing VASP on those input files, and (3) process the resulting VASP output. Also included is a Java™ (Ref. 50) program (ProjectWave.java) that is called by the shell scripts generated by the Mathematica notebook; ProjectWave.java serves to modify the occupations of KS orbitals in WAVECAR files. Details of how to use this code are provided in example form in Sec. IV of the GenMODiagrams.nb Mathematica notebook.

Since the electrostatic force component on each ion due to a given KS orbital $F x i ε$ (defined above in Sec. II) is not a default output of VASP, we calculate it by setting the occupation to zero for all orbitals that do not have an eigenvalue within some small (of order the mean eigenvalue spacing) range of ε in the WAVECAR file (using ProjectWave.java), and by then allowing VASP to calculate the force on each ion due to the remaining charge density (while not allowing VASP to perform electronic steps that would change the charge density). The result of such a calculation necessarily includes the forces imparted by the net positive charge (sum of positive nuclear charge and negative core charge) of the ions. Because we are concerned only with the KS states that we selected to be occupied, $F x i ε$ should include only the force components due to the valence electrons in the calculation and not the contributions from the frozen ionic cores. To address this unwanted counting of ion induced forces, we additionally let VASP calculate the forces for the equivalent system, but with all of the KS orbitals unoccupied, and then subtract off this background. The forces from the fully unoccupied system can then be added back, but with a weight of N_ε/N, where N is the number of occupied electrons of the original system and N_ε is the number of electrons included in the small interval about ε, in order to determine the effective force on each ion, $F ̃ x i ε$⁠.

These values are then used in Eq. (B4), where we make the approximation: $F ̃ x i e , ε =〈ε|e〉〈e|ε〉 F ̃ x i ε ≈〈ε⋮e〉〈e⋮ε〉 F ̃ x i ε$⁠. As a concrete example, 〈e⋮ε〉 would be $〈 ψ m AQ ⋮ ψ n A Q − adsorbed 〉$ in the case of the AQ molecule, where the index m ranges over the KS states of AQ in vacuum and the index n ranges over the near continuum of KS states of AQ adsorbed on a Cu surface. The integral from A to B of Eq. (B4) is then calculated numerically (using the trapezoid rule) over the discrete set of configurations shown in Fig. 8, where A becomes the (a) configuration (hcp) and B becomes the (i) configuration (fcc), yielding as a final result the values of $Δ E ̃ A → B e$ for each MO of AQ, presented in bar graph form in Fig. 6 of the main text. The implementation of the integration scheme for this example can therefore be expressed as

where the integral of Eq. (B4) has been numerically approximated by a discrete sum over interpolation indexes j from 1 to 8 (representing (b) through (i) of Fig. 8); note that the j index on each nuclear coordinate component must be included since the nuclear coordinates vary over each interpolated configuration. The factor of 1/2 and the inclusion of forces at two interpolation configurations in each summand amount to approximating the area under the curve as a trapezoid. In this example, A and B would then be the interpolations (a) and (b), respectively, while e and ε would be $ψ m AQ$ and $ψ n , j A Q − adsorbed$⁠, respectively; the explicit inclusion of the index j emphasizes that the adsorption system orbitals change in each different interpolation configuration.

All of the code along with a subset of the data sufficient to execute the code can be found in Ref. 48.