Design Principles for Metastable Standing Molecules.

Molecular nanofabrication with a scanning probe microscope (SPM) is a promising route toward the prototyping of metastable functional molecular structures and devices which do not form spontaneously. The aspect of mechanical stability is crucial for such structures, especially if they extend into the third dimension vertical to the surface. A prominent example is freestanding molecules fabricated on a metal which can function as field emitters or electric field sensors. Improving the stability of such molecular configurations is an optimization task involving many degrees of freedom and therefore best tackled by computational nanostructure design. Here, we use density functional theory to study 3,4,9,10-perylene-tetracarboxylic dianhydride (PTCDA) standing on the Ag(111) surface as well as on the tip of a scanning probe microscope. We cast our results into a simple set of design principles for such metastable structures, the validity of which we subsequently demonstrate in two computational case studies. Our work proves the capabilities of computational nanostructure design in the field of metastable molecular structures and offers the intuition needed to fabricate new devices without tedious trial and error.

Introduction

The scanning probe microscope (SPM) is a powerful tool for nanofabrication, offering a promising long-term perspective for the prototyping of nanometer-scale functional molecular structures well beyond those accessible by self-assembly.[1−14] Metastable molecular configurations created by SPM-based manipulation can possess unique electronic and mechanical properties, making them interesting for new types of devices. Tip-attached standing molecules in scanning quantum dot microscopy (SQDM)[9,15] are, for example, one of the few actually useful single-molecule devices. Such examples show that a weak coupling between a molecule and its support can, in fact, be a necessary prerequisite for functionality[9,11,15,16] and thus an integral aspect of future molecular devices. When accepting this fact, stability becomes a crucial design issue because weak coupling often implies metastability. Typically, metastable molecular structures have to be fabricated and investigated at cryogenic temperatures of 5 K to suppress instability from thermal fluctuations.[17] If they could, however, reliably withstand higher temperatures, strong electric fields and other effects such as charging, the application range of molecular devices would greatly increase. Here, we demonstrate how computational nanostructure design using density functional theory (DFT) calculations allows us to discover and systematically improve suitable metastable molecular devices.

DFT is the workhorse of contemporary atomistic quantum mechanical calculations and an important complement to experimental observations. Notably, it has provided insights into many metal–organic hybrid systems like molecular monolayers adsorbed on high-symmetry surfaces.[18−26] The capability of DFT to describe such systems has steadily improved, particularly due to a more sophisticated treatment of long-range van der Waals (vdW) interactions.[24,27−31] However, the metastable conformations encountered in SPM-based nanofabrication (Figure ) are substantially different from the monolayer systems that DFT is typically benchmarked against.

Figure 1

Fabrication of standing PTCDA. To bring PTCDA into a vertical conformation on the SPM tip (right) or the adatom surface (left), first a bond between the SPM tip and a Ocarb atom of the surface-adsorbed molecule needs to be established. Upon tip retraction (center), either the molecule–surface connection breaks (right), or if a pedestal of two adatoms has been fabricated at the lower end of the molecule, the tip–molecule bond eventually breaks and the molecule remains vertically on the surface (left).

In a previous paper, we could, nevertheless, prove experimentally that the Perdew–Burke–Ernzerhof (PBE) functional[32] in combination with the many-body dispersion approach (MBD)[28−30] can accurately reproduce the potential energy landscape of a standing molecule.[17] Specifically, the computed barrier height of approximately 30 meV which prevents the collapse of a standing 3,4,9,10-perylene-tetracarboxylic dianhydride (PTCDA) molecule on adatoms on Ag(111) was confirmed by thermal excitation measurements. Since barrier heights are practically energy differences and not absolute energies, the accuracy of the respective calculation is well below typical uncertainties of absolute DFT energies of ∼100 meV. While the 1–2 meV difference between theory and experiment in ref (17) is probably a fortuitous agreement, also larger discrepancies would have no negative effect on the conclusions and design rules derived here. This confirmation of our DFT method for one specific metastable configuration is a crucial prerequisite for the present in-depth study.

The standing PTCDA on Ag(111) is characterized by a competition between stabilizing covalent and destabilizing vdW interactions that causes a shallow potential energy minimum and thermally induced collapse between 10 and 14 K.[17] This insight is our starting point for a deeper and more detailed study culminating in the development of practically applicable design rules. To identify global trends, we examine a broad set of possible configurations of standing PTCDA (s-PTCDA), including cases in which PTCDA is attached to the SPM tip as is required for SQDM.[9,15,16] Due to the variable and experimentally inaccessible atomic structure of the SPM tip, it is not possible to determine the configuration of the PTCDA-tip bonding region unambiguously. By computationally exploring the potential energy surface of s-PTCDA on various tip apex structures, we can nevertheless reveal fundamental aspects of the PTCDA-tip attachment, like the number of involved Ag–O bonds.

Because the potential energy surface of s-PTCDA is the consequence of a competition between covalent and vdW interactions,[17] the design criteria for enhanced stability can be formulated with regard to these two interaction types. We have organized the paper accordingly: After summarizing the computational methods, discussing the experimentally observed properties of s-PTCDA and defining the set of SPM tip structures used in our analysis, we first analyze the chemical bonding of s-PTCDA to the Ag structures and subsequently the influence of the vdW interactions on stability. The deduced design criteria are finally validated in two computational case studies.

Computational Methods

DFT Calculations

We performed DFT calculations using the Fritz-Haber-Institute ab initio molecular simulation package (FHI-aims).[33] We used the semilocal exchange-correlation functional of Perdew, Burke, and Ernzerhof (GGA-PBE)[32] to treat the electronic exchange and correlation. To account for long-range dispersion interactions, we used the vdWsurf scheme[24,27,31] and the MBD method.[28−30] We use vdWsurf when exploring tip structures and MBD for our computational case studies where highly accurate vdW energies are beneficial. To include relativistic effects, we applied the scaled zeroth-order regular approximation (ZORA)[34] in all calculations.

Depending on whether the target system is an isolated cluster tip or a periodic slab, we treated them with either aperiodic or periodic boundary conditions, respectively. For the periodic slab calculations, we used a Γ-centered 2 × 2 × 1 Monkhorst–Pack grid[35] to sample the Brillouin zone during structural relaxation. It was replaced with a denser 4 × 4 × 1 mesh for PBE + vdWsurf or PBE + MBD production calculations and charge density difference analysis. For the MBD postprocessing calculations, the Brillouin zone was sampled with a half-Γ-shifted 8 × 8 × 1 K-point grid. To aid the convergence, we applied a 0.02 eV broadening to all states, using a Gaussian occupation smearing scheme.[36]

Using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm, we relaxed the structures via two serial steps, first by expanding the Kohn–Sham wave function with the default numerical “light” basis sets,[37] applied to all atomic species. Once finished, all basis sets were replaced with the default numerical “tight” settings,[37] and the relaxation was resumed. We continued the relaxation until the maximum force on each atom, in either setup, was less than 10–3 eV/Å. To obtain a well-converged electronic description of the systems, a threshold of 10–7 eV for the total energy, 10–4 eV for the sum of eigenvalues, and 10–6 e/Å3 for the charge density was applied during all SCF cycles.

Tip Model Preparation

To prepare our tip models, we first built a ten-layer-thick 10 × 10 Ag(111) slab with the experimental lattice constant and the surface normal pointing along the z-axis. Subsequently, cluster tips were carved out of the slab by removing silver atoms until the desired shapes were achieved. For the two-atom apex tips, silver adatoms were considered with either the minimum lateral distance in xy, i.e., 2.88 Å (D1) or slightly larger at 3.32 Å (D2). Depending on the atom removal strategy applied to the upper layers, the created clusters resembled either rather sharp or blunt tips. Sharp tips are made by keeping always only the nearest-neighbor (NN) atoms with respect to the previous layer. To create blunt tips, we kept more than just NN atoms in the first layer beneath the apex, depending on how blunt and in what shape the tip model was supposed to be. From the second layer onward, we used the original protocol of keeping just the NN atoms. The single apex Tip 1 is created with six atomic layers (including apex), while all other tips have five layers. Except for the apex atom and the layer above it, all silver atoms in the tip models are frozen during geometry relaxations.

Slab Preparation

To model the Ag(111) surface, we built a four-layer 8 × 8 silver slab consisting of 256 Ag atoms and a 60 Å thick vacuum layer added between slab images along the z direction. Starting from the experimental surface lattice constant (a = b = 2.88 Å), the lattice parameter was first converged for a four-layer Ag(111) primitive unit cell, and the supercell was constructed accordingly. The converged lattice parameter (a = b = 2.875 Å) was very close to the experimental value. The adatom pedestal in the D1 configuration was created by placing two Ag adatoms at neighboring hcp hollow sites 2.36 Å above the surface. Similar to the two-atom apex tips, we defined the reaction coordinate ξ for the rotation around the axis of adatoms as the angle between the plane of PTCDA and the surface normal. The bottom three layers of the slab were held fixed, while the top layer and the adatoms were allowed to relax freely.

One-Bond Scenario Configuration Mapping

In the one-bond scenario, the (azimuthal) angle α is measured between the short axis of PTCDA and the x axis, and the angle β is measured between the long axis of PTCDA and the z axis (Figure b). We map the potential energy surface for both tip models on a grid of 3 α values (0°, 90°, and 180°) and 21 β values in the interval β = [−90°...90°], with (α = β = 0°) denoting the strictly vertical alignment of PTCDA shown in Figure b. Configurations for which no stable constrained relaxation was possible, e.g., due to molecule–metal overlap, were, however, not considered further.

Figure 2

Tip models. (a) All tip models used in this study. The tip apex atoms which bind to PTCDA are colored black for visual enhancement. The flat area of each tip onto which the apex atoms are placed is highlighted. Tip 1 and Tip 2 feature a single apex atom, and all other tip models feature two apex atoms. (b) The geometries of PTCDA attached to a single-apex tip are defined by the angles α and β. Shown is the case of α = β = 0 in a space-filling model of Tip 1 (right) and in a schematic view along the positive z axis (left). (c) In the two-bond scenario, the geometries of PTCDA are defined by the angle ξ. Shown are exemplary geometries for PTCDA on Tip 3 (left), bare Ag(111) (center), and the adatom surface (right). The tip models 3–10 in panel (a) are oriented such that positive ξ values always correspond to a tilting of the molecule to the right.

Constrained Optimization

For a proper mapping of the potential energy surface, structural relaxation is required while keeping the angles α and β of the molecule fixed. In the two-bond scenario, only ξ needs to be considered, and the constraint was imposed by fixing the z coordinate of a single carbon atom at the side of the PTCDA molecule far from the Ag–Ocarb bond. Since the relaxation of the Ag–Ocarb bond leads to a small change in ξ under this constraint, the value of ξ was remeasured once the structure was optimized, and the measured value was used in all plots. In the one-bond scenario, fixing α and β angles was achieved by constraining one coordinate of each of three carbon atoms of the perylene core at the side of the PTCDA molecule opposite of the Ag–O bond. Again, after optimization, the tilt angle β was remeasured.

Results and Discussion

Fabrication and Properties of Standing PTCDA

The design principles formulated in this theoretical study are validated by their capability to predict and explain experimentally observed properties of s-PTCDA. Particularly, the different behavior of s-PTCDA either on the Ag(111) surface or on the Ag-covered SPM tip (Figure ) is of relevance in this respect. Here, we briefly describe the SPM-based fabrication process of each configuration and their experimentally observed properties.

Standing PTCDA on Ag(111) can be created by two-contact manipulation with the SPM tip,[11] where one contact is between the molecule and the tip and the second contact between the molecule and the surface. First, two Ag adatoms which later form a “pedestal” are attached to the two carboxylic oxygen atoms (Ocarb) at the short side of PTCDA by lateral manipulation. Subsequently, PTCDA is contacted by the tip at one of the other two Ocarb atoms and lifted into a vertical orientation. Further tip retraction breaks the tip–PTCDA bond, finalizing the assembly process (Figure , left). If, on the other hand, no Ag adatom pedestal is used and lifting of an isolated PTCDA molecule is performed (Figure , center), the molecule is detached from the surface and typically remains on the SPM tip in a likewise vertical orientation (Figure , right),[15,16,38−40] the detailed atomic configuration of which is experimentally inaccessible.

With this fabrication process, s-PTCDA is a prototypical SPM-fabricated metastable molecular configuration with properties that make it also a single-molecule device: On the Ag(111) surface, it functions as a single electron field emitter[11] and GHz oscillator.[17] Tip-attached s-PTCDA, on the other hand, enables a new microscopy technique, SQDM.[9,15] All these functionalities are enabled by the metastable vertical configuration of s-PTCDA, which is characterized by a minimal overlap between the molecular π-system and tip or sample. Due to this minimal overlap, certain electronic states of the molecule remain decoupled from the metal and can be gated by applying a voltage to the tip–molecule–surface junction. At sufficiently high gate voltage, s-PTCDA can then acquire an additional positive or negative elementary charge[9,11,15,16] which is the fundamental mechanism enabling both SQDM and field emission from s-PTCDA. The weak mechanical coupling between s-PTCDA and the metal, on the other hand, allows the oscillation of the entire molecule with a remarkably low GHz frequency.[17]

The stability-related phenomenology of s-PTCDA on the tip and the surface becomes richer due to the option of gating: While negatively charged molecules remain stable, positive charging tends to destabilize the vertical molecular orientation even at our base temperature of 5 K. Attempts to charge a tip-attached molecule positively occasionally result in a collapse of the molecule onto the tip, whereas surface-standing PTCDA always collapses when positive charging is attempted.[11] These observations suggest that positive charging weakens the stabilization and sometimes completely eliminates the stabilizing local potential energy minimum. The molecule’s variable behavior in response to positive charging moreover implies a dependency of the stabilization potential on the atomic configuration of the tip apex since this configuration is the only property that varies between individual experiments. We thus hypothesize that the stabilization potential of s-PTCDA on the Ag tip is typically larger than that on the adatom pedestal on Ag(111) (called “adatom surface” from now on) on which molecules always collapse upon positive charging. Finally, experiments attempting the creation of s-PTCDA on the bare Ag(111) surface turn out to be unsuccessful.

To explore this broad phenomenology computationally, we have to capture the corresponding variations in the Ag structure of both the tip and the surface. Especially the unknown tip geometry presents a formidable challenge which we address in the following section.

Tip–Molecule Bonding Geometries

Since there is only little experimental information about the SPM tip apex configuration, the actual bonding mode of PTCDA to the tip is unknown. On the most elementary level, this uncertainty boils down to the question of how many simultaneous Ag–O bonds are required to stabilize s-PTCDA on the tip. Two scenarios seem plausible: First, a metastable conformation featuring only a single chemical bond (Figure , right) as suggested by the lifting process in which the surface-adsorbed PTCDA is contacted with the SPM tip at one of its Ocarb atoms and lifted off the surface. Second, a molecule attachment with two Ag–O bonds involving two neighboring Ocarb atoms, similar to the known vertical conformation of the molecule on the adatom pedestal on Ag(111), which features two Ag–Ocarb bonds.[11,17] (Figure , left).

To investigate the two scenarios in detail, we examine the bonding of PTCDA to a variety of different tip models (Figure ). This allows us to determine, for example, how much the immediate surrounding of the tip–molecule binding site influences the bond properties. To study the one-bond scenario, we designed two pyramidal single-atom apex tip models (Tip 1 and Tip 2). For the two-bond scenario, we designed a set of eight tip models (Tip 3–Tip 10), all of which feature a flat Ag(111) area containing 5 ≤ N ≤ 40 atoms (highlighted in Figure a) onto which two adatoms are placed as equivalent tip apexes. Since these models are effectively differently sized and shaped cutouts of an adatom–surface slab, their properties naturally converge to the adatom–surface case for very wide tips (e.g., Tip 10 in Figure a).

Importantly, the two-atom apex tip models, Tip 3–Tip 10, are not necessarily “double tips” in the meaning of the term used in SPM, where it describes tips with two apex atoms at the same distance from the surface. Our tip models would only operate as double tips in SPM if their flat areas were exactly parallel to the imaged surface. This is not likely since, in experiment, there would be no preference for any specific crystal orientation of the tip apex region. Hence, our tip models are not in conflict with the experiments in which the use of double tips was avoided.

The bonding scenario also determines the degrees of freedom of the tip-attached molecule since rotation is possible around three axes in the one-bond scenario (Figure b), whereas there is only a single rotation axis in the two-bond scenario (Figure c). Mapping of the potential energy landscape thus requires two different approaches which are outlined below, starting with the more complicated one-bond scenario.

The configuration space of s-PTCDA attached to a single apex tip is primarily spanned by rotations of PTCDA around the apex atom, which we describe in a coordinate system that is fixed to the tip models and aligned with their high-symmetry directions (Figure b). The rotational and mirror symmetries of Tip 1 and Tip 2 allow limiting the conformation mapping to 1/6 of the half sphere below the apex. A rotation around the z-axis changes the molecule’s azimuthal orientation, while rotations around the x or y axis bring it closer to the tip. To reduce the number of computed configurations, we map this “upward” rotation only in the high-symmetry x–z plane (dark pink in Figure b), effectively eliminating one rotational degree of freedom (around the x axis). The mapping thus proceeds by rotating PTCDA first around the z axis (angle α) and subsequently around the y axis (angle β) for three α values (0°, 90°, and 180°) and a series of β values with Δβ = 10°. Small translations that optimize the geometry of the Ag–Ocarb bond are accounted for by constrained structural relaxation (see Computational Methods). Since our constraining method allows also for small relaxations in β, the angle is remeasured after relaxation.

In the two-bond scenario, two adjacent Ocarb atoms of PTCDA are bound to the two Ag apex atoms, which reduces the molecular degrees of freedom to a rotation around the respective Ocarb–Ocarb axis (tilt angle ξ) and small translations accounted for by constrained geometry relaxation (Figure c). Since our two-apex tip models are cutouts of a Ag(111) slab, we place the apexes as adatoms in the known D1 or D2 configurations.[11,17] In D1, the atoms are located in two adjacent hollow sites of the same type (fcc or hcp sites, respectively), while D2 features adatoms in nonidentical hollow sites (fcc–hcp combination). Previously, we found that the D2 configuration provides slightly higher stability to s-PTCDA on the Ag(111) surface, specifically a 5 meV higher potential energy barrier against collapse.[17] Here, six of our tips feature a D1 adatom configuration, while Tip 7 and Tip 9 feature a D2 configuration.

To put the results of our calculations for the various double-apex tip models into perspective, we also include the case of s-PTCDA on the adatom surface for which the structure is known and the hypothetical case of s-PTCDA standing on the bare Ag(111) surface (Figure c). Both configurations feature a single rotation axis similar to the two-apex tips. As motivated earlier, the first step of our analysis will now be a discussion of the (typically stabilizing) DFT-PBE potential energy surface, followed later by an analysis of the (destabilizing) long-range vdW potential.

Contributions of Covalent and Electrostatic Interactions to Molecular Stability

It was recently shown that the covalent and electrostatic interactions captured by DFT calculations with the PBE functional[32] stabilize s-PTCDA on the adatom surface against collapse.[17] In the following, we extend this picture by mapping the PBE potential energy landscape EPBE of s-PTCDA attached to the variety of tip models described above. Our goal is to identify the most likely attachment scheme of PTCDA to the tip after loss of contact to the sample surface and to identify universal design principles that allow optimizing the stability of the vertical state. To obtain EPBE = EPBEsys – EPBEtip – EPBEmol, we subtract the energy of the isolated systems of the tip (EPBEtip) and molecule (EPBEmol) from the energy of the combined system consisting of the tip and molecule (EPBEsys). Although we do not yet discuss the long-range vdW energies EvdW at this point, vdW interactions (calculated with the vdWsurf method) have, of course, been considered in the constrained geometry optimizations for all configurations presented below.

For each tip model, we plot the relative PBE energy ϵPBE(α,β) = EPBE(α,β) – EPBE(0°, 0°) or ϵPBE(ξ) = EPBE(ξ) – EPBE(0◦), respectively, in Figure a and b. It shows that there is a clear qualitative difference between the one-bond and the two-bond scenarios. In the two-bond scenario in Figure b, the stabilizing potentials are quite similar across all tip and surface models, featuring a minimum at ξ = 0° and torsional spring constants between κ = 3 × 10–20 Nm and κ = 5 × 10–20 Nm. This clearly demonstrates that the stabilizing potential must be primarily determined by the two Ag–Ocarb bonds and their immediate vicinity, which is the same for all two-apex tip models. The overall tip shape, on the other hand, seems to cause rather weak modifications of the potential, which could arise, for example, from electrostatic interactions between partial charges on PTCDA and Smoluchowski dipoles[41] of undercoordinated metal atoms at nearby edges and corners of the tip models.

Figure 3

DFT-PBE potential surface of s-PTCDA. (a) PBE energies for the single-atom apex tip models (Tip 1 and Tip 2) as a function of azimuthal molecule orientation α and tilt angle β. The geometries corresponding to the lowest and highest β value are displayed in panel (c). Spline fits of the data are a guide to the eye. (b) Potential energy profiles for s-PTCDA on all double-apex tip models, the bare and the adatom surface as a function of tilt angle ξ. For the corresponding geometries, see Figure a,c. (c) Geometries with minimal and maximal β values for the six curves in panel (a).

Surprisingly, our results in Figure b also show that s-PTCDA experiences a stabilizing PBE potential even when placed on the bare Ag(111) surface. This result suggests that no undercoordination of the metal atoms participating in the Ag–Ocarb bonds is required for stability. However, undercoordination still plays a central role in the actual manufacturing of s-PTCDA by SPM. The high strength of the Ocarb–adatom bonds with a PBE energy of EPBE = −1.05 eV allows breaking the bond between the tip and s-PTCDA by simple tip retraction (Figure left), whereas the much weaker bonds between Ocarb and Ag atoms in the surface (EPBE = −0.29 eV) break upon tip retraction, leaving the molecule on the tip (Figure , right). The configuration of a vertical PTCDA on bare Ag(111) and without contact to the tip would therefore be hardly accessible in an SPM manipulation experiment because it would have to involve breaking the tip–molecule bond, for example, by inelastic excitation of molecular vibrations via a tunneling current.

While all configurations in the two-bond scenario provide a similar level of stabilization, the situation is less clear in the one-bond scenario in Figure a. With only a single molecule–metal bond, there are many nonequivalent trajectories along which the molecule can tilt toward the tip. First, our calculations reveal that the vertical molecule (β = 0°) can rotate around the z-axis practically without encountering a potential energy barrier since, at β = 0°, the energies for α = 0°, α = 90°, and α = 180° differ by less than 10 meV. Furthermore, it shows that the PBE potential ϵPBE(α,β) does not generally stabilize the vertical molecular orientation, as it does in the two-bond scenario. Instead, the potential energy curves of Tip 2 for α = 0◦ and α = 180◦ have minima around β = +25° or β = −30°, respectively, and the β = 0° minimum for Tip 1 at the α = 0◦ orientation is negligibly shallow. Other tilt paths, on the other hand, do provide restoring forces with path-dependent torsional spring constants between κ = 2.6 × 10–20 Nm and κ = 5 × 10–20 Nm, similar to the two-bond scenario.

The variability of the potential energy profiles in the one-bond scenario is caused by the tendency of the adjacent (initially unbound) Ocarb atom to form a Ag–O bond as well. This transforms the one-bond scenario into a two-bond scenario as depicted in Figure c (α = 180° for β < 0 and α = 0° for β > 0). Since the energy scale of bond formation substantially exceeds the energy scale of mere bond tilting, bond formation strongly deforms the potential energy surface and is typically always favored. The energy gain resulting from the second bond depends significantly on the Ag environment to which the second Ocarb atom binds. This is exemplified by the EPBE values for the bare and the adatom surface quoted above and by the different depths of the minima at (0◦, 25◦) and (180◦, – 30◦) of Tip 2. As in the two-bond scenario, the remaining variations in the ϵPBE curves are likely related to inhomogeneous electrostatic potentials from Smoluchowski dipoles. For Tip 1, the second bond can only be established at rather large tilt angles β. At these geometries, our approach of constraining only a single C atom in the x direction and two C atoms in the y direction to fix both β and α while allowing a maximal freedom of structural relaxation (see Methods) fails insofar as the molecule adopts configurations where it does not bind to the apex atom anymore but favors binding to other atoms at the edge of the tip. In terms of a rigorous analysis, we have excluded the respective geometries, which does not affect any of our conclusions.

We obtain two central results from our calculations of the PBE energy EPBE of s-PTCDA attached with one or two Ocarb atoms to SPM tip models of various shapes. On one hand, we find that Ag–Ocarb bonds indeed stabilize the vertical orientation of s-PTCDA such that, in the hypothetical absence of all vdW interactions, the two-bond scenario would generally be stable. The single-bond scenario, on the other hand, is a very unlikely outcome of a PTCDA lifting experiment since there will usually be the option to form an energetically favorable second Ag–Ocarb bond. The possibility of the molecule to rotate freely around α at β = 0° is important in this respect since it provides access to pathways along which this second bond can form without the necessity to overcome a potential energy barrier first. Hence, while the lifting of PTCDA always starts with a single Ag–Ocarb bond, our calculations suggest that a second bond to the tip will practically always form at some point of the lifting process, at the latest when the molecule is fully detached from the surface.

Point Charge Model of the Stabilization Mechanism

While the PBE energies clearly show the presence of a stabilizing potential, a deeper analysis is required to rationalize this finding and to formulate a simple design rule for increased stability. When s-PTCDA is placed on the adatom pedestal, the bonding causes a charge redistribution which can be visualized in a charge density difference plot (Figure a). Interestingly, this redistribution further increases dipoles which already exist in the isolated systems, namely, the polar C–O bonds of PTCDA and the Smoluchowski dipoles of the Ag adatoms (Figure a, right).[15] In the case of s-PTCDA on the bare Ag(111) surface, dipoles are created in the Ag–Ocarb bonding region as well (Figure a, left), despite the absence of dipoles in the isolated Ag(111) slab. Thus, both cases, with and without an adatom pedestal, are characterized by dipoles which point away from the surface. Importantly, the adatom Smoluchowski dipoles and even the Ag–Ocarb bond dipoles on the bare surface do not change their orientation when the molecule is tilted, while the C–O bond dipoles obviously do (Figure a).

Figure 4

Electrostatic model of the stabilization mechanism. (a) Charge density difference Δn(r) = nsystem(r) – [nmolecule(r) + nslab+adatoms(r)] for PTCDA on bare Ag(111) (left) and on the adatom surface (right). The isosurfaces are displayed for −0.005 eÅ–3 ≤ Δn(r) ≤ +0.005 eÅ–3 and negative charge flows from red to blue regions. Local dipoles which are created or enhanced by the charge redistribution are indicated by yellow arrows. (b) Two-dipole model for the example of PTCDA on the adatom surface. Partial charges from Smoluchowski dipoles and bond dipoles are indicated in the space-filling model. Charges colored in gray represent the Ag adatom dipole; red charges represent the O atoms; and turquoise charges represent the nonaromatic C atoms. The angle ξ′ equals the PTCDA tilt angle. (c) Lowest-energy configurations of the two dipoles for a series of tilt angles ξ′. (d) Properties of the two-dipole model for r1 = r2 = d = 1 Å and q1 = q2 = 1e. Red and blue symbols show the optimal (i.e., lowest energy) θ(ξ′) relation for the dipole model (blue squares) and for the DFT calculation (red circles). The lines show the (relative) energy of the respective configurations (dipole model: blue, DFT: red). A scaling factor of 1/2 has been applied to the DFT data (see text). The inset shows the Coulomb energy Eϕ(θ, ξ′) (eq ) where regions with Eϕ > 0 are colored black. (e) Lowest-energy potentials ϵϕ(ξ′) for a net charge δq which offsets each of the two charges that form the upper dipole pC–O. δq = 0 reproduces the blue curve from panel (d).

Since we have already established that the stabilization must originate from the immediate surrounding of the Ag–Ocarb bonds, observing such local charge (re)distributions raises the question whether and how they are involved in the stabilization mechanism. To study this aspect, we set up a simple analytical model in which the Ag-PTCDA bonding region is solely described by the interaction between four charges forming two dipoles (Figure b). The dipole which mimics the polar Ag–Ocarb bonds (pAg) is fixed vertically to the surface, while the second one mimicking the polar C–O bonds (pC–O) has a variable orientation. More rigorously, the second dipole represents the projection of the C–O bond dipole onto the long molecular axis since all four point charges in our model are confined to one plane (Figure b).

The model is parametrized by seven quantities (r1, r2, d, θ, ξ′, q1 > 0, and q2 > 0), and the electrostatic potential energy of this charge arrangement can be expressed as

For our semiquantitative analysis, we use a parametrization with generic yet plausible values of r1 = r2 = d = 1 Å and q1 = q2 = 1e. The angles θ and ξ′ have direct equivalents in the DFT calculations since θ corresponds to the angle between the Ag–Ocarb bond and a plane that is normal to the surface and passes through both Ag adatoms, while the tilt angle ξ′ of the upper dipole pC–O represents the molecule tilt angle ξ. The local lateral relaxation of the Ag–Ocarb bonds at a given ξ which we perform in DFT corresponds to choosing an angle θ which minimizes the Coulomb energy of the two-dipole system (Figure c). The respective function Eϕ(θ, ξ′) (inset of Figure d) has minimal values along an almost linear θ(ξ′) relation, which thus marks the minimal-energy path along which the Ocarb atom moves in the two-dipole model upon tilting of the molecule (pC–O dipole).

Two properties of the two-dipole model can be compared to our DFT results, namely, the outcome of the described relaxation expressed in the relation θ(ξ′) and the ξ′-dependent relative Coulomb energy ϵϕ(θ(ξ′), ξ′) = Eϕ(θ(ξ′), ξ′) – Eϕ(0, 0) of the two-dipole system. The calculated θ(ξ′) relation very well matches the corresponding relation of the full DFT calculation in the relevant range of −60° < ξ < 60° (blue and red symbols in Figure d). Since, in the model, this relation does not depend on the values of q1 and q2 at all and is moreover rather insensitive to the choice of r1, r2, and d, the good correspondence to the DFT data is a robust conclusion. This indicates that the electrostatic forces as described by the two-dipole model could (at least partially) be the stabilization mechanism of s-PTCDA.

The plausibility of this assumption is further underpinned by a comparison of the tilt-angle-dependent relative Coulomb energy ϵϕ(ξ′) (Figure d, blue curve) to the relative PBE-DFT total energy ϵPBE(ξ) (Figure d, red curve). Here, we have scaled ϵPBE(ξ) down to match the ϵϕ(ξ′) curve in the range of ±30° and found that a factor of 1/2 was required. This is a reasonable result, given that there are two adatoms and Ag–Ocarb bonds and two C–O bonds in the s-PTCDA DFT calculation but only one dipole for each in the two-dipoles model. This rather good correspondence between ϵϕ(ξ′) and ϵPBE(ξ) across a wide range of tilt angles also implies that our choice of model parameters r and q was reasonable and that the two-dipole model can predict the correct order of magnitude of the stabilizing potential.

Finally, we turn to the experimental observation that a positive charging of s-PTCDA by appropriate gating destabilizes its vertical configuration,[11] particularly in the case of s-PTCDA on the adatom surface, where the stabilizing potential is rather shallow (Figure d). A negative charging of the molecule, on the other hand, does not cause any destabilization.[9,15,16] Using the two-dipole model, this observation can be attributed to an electrostatic effect as well. To mimic the charging, a small net charge δq is added simultaneously to both of the charges which form the dipole pC–O that represents the molecule, such that they become +q2+δq, and −q2+δq. As a consequence, the ϵϕ(ξ′) curve is flattened and finally inverted for δq > 0 (i.e., positive charging), while it becomes steeper for δq < 0 (negative charging) (Figure e), in correspondence to the experimental observations. This behavior emerges because tilting decreases the distance of the molecule charges +q2+δq and −q2+δq from −q1. If δq is negative, this increases Eϕ, while positive δq decreases Eϕ.

Figure 5

Relative vdW and total potentials in the one-bond and two-bond scenario. (a) vdWsurf potentials ϵvdW(β) in the one-bond scenario. (b) vdWsurf potentials ϵvdW(ξ) for the two-bond scenario. (c) Relative PBE + vdWsurf total energy ϵ(β) for the one-bond scenario. (d) Relative PBE + vdWsurf total energy ϵ(ξ) for the two-bond scenario. Spline fits of the data are a guide to the eye.

Because the two-dipole model is based on point charges, it does not, for example, account for mutual (de)polarization effects of the dipoles. Nevertheless, given its simplicity, the correspondence to the DFT calculations and experiments is convincing. In essence, the tilting of the upper one of two stacked and initially collinear dipoles provides the restoring force which stabilizes the vertical molecular orientation. With this explanation, we can derive a simple design rule aimed at engineering the stabilizing potential: Since at least a part of the stabilization emerges from similarly oriented dipoles in the molecule and the metal, the selective modification of these dipoles is the key to a modification of stability. We will show an example for such an engineering of electrostatics when demonstrating the application of our design rules in a later section.

Destabilizing van der Waals Interactions

Here, we turn to the long-range vdW interaction which provides the second contribution to the total potential energy of the molecule–metal structures and ultimately decides whether a configuration considered stable according to its PBE potential will indeed be stable. The vdW attraction between the molecule and the surface was found to cancel the stabilizing PBE potential for s-PTCDA on the adatom surface almost completely, allowing a thermally induced collapse already at T = 10 K.[17] In the following, we discuss the general impact of vdW forces on the stability of standing molecules using our tip models and the bare and the adatom surface (Figure a,b). As in the case of covalent interactions, only relative vdW energies ϵvdW = EvdW – EvdW(0°) are relevant for the stability of s-PTCDA, where EvdW(0°) is the energy of the vertical configuration.

Depending on the context and aim, we either use the vdWsurf method[24,31] based on screened pairwise interactions or the MBD method,[28,29] which implicates a more accurate nonlocal description of vdW interactions at a higher computational cost. We use vdWsurf when exploring tip structures, an application in which vdW interactions are primarily determined by the atomic configuration of the tips and not by the difference between the two computational methods. We use MBD for our computational case studies (below) where the structures are precisely known and highly accurate vdW energies are beneficial.

Since the long-range metal–molecule vdW interaction is generally attractive, it does not provide any stabilization to s-PTCDA. In that sense, the data in Figures a and b can be interpreted rather straightforwardly in terms of changes in the distance between the tip or surface and PTCDA. In the two-bond scenario, the potential energy maximum is always at ξ = 0° because any tilting of the molecule will bring all its atoms closer to the tip cluster or the surface (Figure b). The situation is different in the one-bond scenario, where tilting at (α = 0°, β < 0°) and (α = 180°, β > 0°) causes some C and O atoms in the vicinity of the Ag–Ocarb bond to initially withdraw from the tip, such that (shallow) potential energy maxima are found at nonzero β values instead (Figure a). Besides the general lack of stabilization, the variations among the individual ϵvdW(β) and ϵvdW(ξ) curves are quite pronounced. Since our computational method vdWsurf is, in essence, a pairwise approach, these variations can be interpreted straightforwardly by the Ag atom count in the tip models and by the tip sharpness. Sharper tips generally cause less vdW interactions as they have fewer atoms in the vicinity of the molecule. Consequently, s-PTCDA on the adatom surface and on the bare Ag(111) surface have the strongest vdW interactions since their Ag atom count is high and they are the analogue of completely blunt tips. Since this aspect is independent of the specific molecular species, it is valid for a broad range of SPM tip functionalizations.

Since we have studied a wide variety of tip models, we can estimate the region of influence around the molecule within which the tip structure has an impact on the stability of s-PTCDA. In Figure , we therefore compare the relative vdW energies for an exemplary tilt angle of ξ = 35° for a range of tip “surface” areas (highlighted regions in Figure a) in the two-bond scenario. It shows that our biggest tip model (Tip 10) with 40 atoms in the “surface” plane already provides almost the same vdW energy gain upon tilting as the infinitely extended adatom surface. Hence, we can conclude that regions of the tip which are further from the apex atom than the length of the standing molecule itself have negligible influence on the stability.

Figure 6

Relative vdWsurf energy ϵvdW(ξ = 35°) as a function of the number N of Ag atoms in the surface layer of the different tip models (highlighted area in Figure a). For the adatom surface, N = ∞.

A crucial result is found when comparing s-PTCDA on the adatom surface and on the bare Ag(111) surface (Figure b): In the absence of the adatom pedestal, the vdW interaction increases much more strongly when tilting the molecule, which is a consequence of the reduced molecule–surface distance. In combination with our results for the PBE potential, we can thus conclude that the primary role of the adatom pedestal is not to enhance the stabilizing part of the potential but rather to act as a simple spacer separating PTCDA from the Ag atoms of the surface. This separation reduces the destabilizing effect of the vdW interaction considerably and allows s-PTCDA on the adatom surface to be stable. Besides its role as a spacer, the pedestal is, however, also essential for the practical fabrication of s-PTCDA on the adatom surface by SPM manipulation (see above).

As expected, the variations in ϵvdW also cause strong variations in the total relative energies ϵ = ϵPBE + ϵvdW among the different tip models (Figures c,d), such that the depth of the stabilizing potential varies by about 1 order of magnitude between Tip 3 and the adatom surface. This is in general agreement with the experimental finding that s-PTCDA on the SPM tip is typically more stable and can often also survive positive charging.[9,15] Only on the bare Ag(111) surface do the destabilizing vdW interactions completely outweigh the stabilizing ϵPBE potential such that, even if it could be fabricated, s-PTCDA would not be stable on bare Ag(111) (Figure d). Moreover, the tilt paths along which the transition from a one-bond to a two-bond scenario occurs for Tip 1 and Tip 2 are clearly visible in the total relative energy plot (Figure c).

From our analysis of the long-range vdW interactions, we have obtained a second design rule for increased stability of metastable molecular structures on metals: The vdW attraction should be reduced as much as possible, which can happen either by increasing the distance between molecule and metal, as in the case of sharp tips or the adatom surface, or by reducing the number of atoms in the molecule. In the final part of our paper we will now close the loop and demonstrate that the deduced design principles indeed allow a systematic engineering of stability.

Computational Case Studies

Here we validate the two basic design principles that we have obtained in this study. Since our calculations for a variety of tip models have already provided ample examples for vdW engineering by modification of the metal structure, we will instead examine a reduction of the molecule size. The engineering of dipoles at the molecule–metal bonding region will, on the other hand, be exemplified by changing the atomic species of the pedestal adatoms.

Consequently, we now examine the relative energy profiles ϵ(ξ) = ϵPBE(ξ) + ϵvdW(ξ) of two configurations which are both adaptations of s-PTCDA on the adatom surface (Figure a). On one hand, we decrease the vdW interaction by replacing PTCDA by the smaller NTCDA molecule which lacks one naphthalene unit but has the same functional groups. On the other hand, we engineer the dipoles in the bonding region by replacing the Ag adatom pedestal by a pedestal of Pt adatoms. Given the high Pt(111) work function of almost 6 eV,[42] Pt adatoms on Ag(111) are expected to create strong dipoles pointing toward the surface, while Ag adatoms on Ag(111) exhibit a rather strong dipole of 0.66 D pointing away from the surface.[15] According to our first design rule, this should destabilize PTCDA on a pedestal of Pt adatoms.

Figure 7

Engineering of the stabilization potential. (a) Structures of s-PTCDA on Ag (top) and Pt (bottom) adatoms and s-NTCDA on Ag adatoms (center) on D2 adatom pedestals. (b) Relative energy profiles ϵPBE (left), ϵvdW using MBD (center) and ϵ = ϵPBE + ϵvdW (right) as a function of tilt angle ξ for the three examples in panel (a). The ξ = 0° values for PTCDA on Pt are only estimated since the respective configuration is unstable under our minimal constraining.

Relative PBE and vdW energies (using the MBD method) for a series of tilt angles ξ are shown in Figure b. Indeed, it shows that the stabilizing nature of the covalent Ag–Ocarb bond (ϵPBE) is barely altered when replacing PTCDA by the chemically similar NTCDA on the adatom surface. Using Pt adatoms beneath PTCDA, however, changes the picture completely. As expected from the two-dipole model, the oppositely oriented dipoles pPt and pC–O now have a destabilizing effect on s-PTCDA, initially leading to a reduction in ϵPBE as the molecule is tilted. The DFT calculation, however, also reveals the limits of our simple point charge model which does not predict the observed increase in ϵPBE beyond |ξ| = 30°. This increase points toward a stabilization mechanism with two energy contributions: While we have turned the Coulomb interactions from stabilizing to destabilizing when switching from Ag to Pt adatoms, a second stability mechanism, possibly related to exchange interactions, is still in place. Irrespective of the precise mechanism, our design choice has practically removed the stabilizing effect of the ϵPBE(ξ) potential as intended, thus validating our first design principle. The presence of two stabilization mechanisms and the corresponding double-well potential (Figure b) could, on the other hand, even be an opportunity allowing the fabrication of molecules in tilted configurations if the vdW interaction is sufficiently small.

The MBD vdW energies in Figure b validate also the second design principle, since replacing PTCDA by NTCDA yields the expected decrease in ϵvdW. The magnitude of this decrease of about 20–25% is, however, considerably below what would be expected from the heavy atom count alone (20 for NTCDA compared to 30 for PTCDA). This fact, again, illustrates the strong distance dependence of the vdW interactions: While NTCDA is smaller, its center of mass is also closer to the Ag(111) surface, partially compensating the size-related effect. The (unintended) reduction of ϵvdW for PTCDA on the Pt adatoms can be explained by the same effect. During tilting, PTCDA on Pt approaches the surface slower than PTCDA on Ag adatoms at the same tilt angles (this would correspond to a flatter θ(ξ) relation in Figure ), which reduces the vdW interaction for PTCDA on Pt.

Ultimately, the stability of both scenarios, NTCDA on Ag and PTCDA on Pt adatoms, can only be judged based on the relative total energies ϵ = ϵPBE + ϵvdW plotted on the right of Figure b. It shows that the strategy derived from our design rules proved successful. NTCDA on the Ag adatom surface has a higher overall stabilization barrier than PTCDA on the same adatom type, while PTCDA on Pt adatoms has no stabilizing barrier at all. While this second case resembles PTCDA on bare Ag(111) (Figure d), the similarity is only superficial: In PTCDA on bare Ag(111), the destabilizing vdW interaction is increased, while for PTCDA on Pt adatoms the stabilization by electrostatic and covalent interactions is eliminated.

Conclusions

SPM-based fabrication is a powerful tool which grants access to unexplored molecular configurations with interesting properties and can draw from an almost infinite pool of molecular building blocks. Experimentally searching this vast chemical and configuration space for structures with a desired functionality can be extremely tedious if no clear strategy is available. Here, we have shown that computational nanostructure design can provide such strategies and straightforward design rules. The results of our study are relevant for experiments and computation alike: The design rules for standing molecules can be directly employed to conceive future experiments, while our demonstration that computational nanostructure design is well suited to describe the potential energy profile of molecular nanostructures may trigger its application to new device ideas beyond standing molecules. Here, specifically the insight into the role played by dipole–dipole interactions in nanostructures is of very general nature. While we have applied our computational methods to the specific topic of metastable standing molecular configurations, it can probably be extended to other potential-energy-related aspects of molecular nanostructures and would also offer clear design rules. Studying the specific topic of metastability has the benefit that this aspect is tightly linked to functionality. This is exemplified by the role of s-PTCDA as a sensor, field emitter, and mechanical oscillator. Moreover, owing to their properties as well-defined, weakly coupled quantum systems, we also envision standing molecules to play a role in future quantum devices, for example, quantum sensors. The finding that PTCDA practically always attaches with two bonds to the SPM tip if retracted sufficiently from the surface not only has deepened our understanding of the role of s-PTCDA as a sensor in scanning quantum dot microscopy but also is important for SPM-based manipulation in general because the tip–molecule relation is that of the actuator and the work piece, and the presence of two bonds reduces the degrees of freedom of the latter. Summarizing, our findings mark an important step toward exploiting the full potentials of nanofabrication with scanning probe microscopes.