Controlling Adsorbate Diffusion on a High-Symmetry Surface through

ARTICLE pubs.acs.org/JPCC

Controlling Adsorbate Diffusion on a High-Symmetry Surface through Molecular Shape Selection David M. Huang* School of Chemistry & Physics, The University of Adelaide, SA 5005, Australia

Peter Harrowell School of Chemistry, University of Sydney, NSW 2006, Australia ABSTRACT: Using a simple rigid model of a chemisorbed molecule, we explore the relationship between a molecule’s shape and its diffusion on a (111) surface of a face-centeredcubic crystal. Anisotropies in the diffusion constant of up to a factor of 15 are found, along with the enhancement of diffusion perpendicular to an applied field and an approximate decoupling of the size of the diffusion barrier from the surface binding energy. The model reproduces experimental measurements of adsorbate diffusion both qualitatively and with reasonable quantitative accuracy yet is simple enough that all possible diffusion pathways can be enumerated with little computational expense.

1. INTRODUCTION A molecule’s shape embodies its fundamental chemical information. Catalysis, molecular locomotion, and signal transduction (e.g., smell) are examples of how aspects of that information can be manifest through the interaction of molecules with a suitably structured substrate. In this paper we consider a simple and accessible example of this principle—the role of molecular shape in determining a molecule's diffusion when chemisorbed on a close-packed crystal surface. The simplicity of the substrate allows us to explore the role played by the adsorbate’s shape and its interaction with the surface. The attachment of molecules to metal surfaces is a staple of nanotechnology fabrication,1,2 and understanding how molecular shape can influence the binding energy and mobility is important. By way of example, a nanotube-electrode junction between a bucky tube attached to a gold surface by thiol groups was found in an ab initio study to exhibit signicant mobility in spite of the five covalent binding sites employed.3 Such mobility can involve marked anisotropy. Kwon et al.4 observed unidirectional motion of 9,10-dithioacetylanthracene (DTA) (see Figure 1) along each of the three high-symmetry [110] axes of a Cu(111) surface. The high symmetry of the substrate here is significant. It is, after all, not surprising that anisotropy in a surface can produce anisotropic (even unidirectional) motion of adsorbate molecules across the surface. This phenomenon has been observed both in theoretical studies of model systems5 and in experiments with highly anisotropic surfaces6
8 such as the (110) surface of fcc metal crystals. What is interesting in the observations of ref 4 is that it suggests that control of adsorbate motion might be exerted r 2011 American Chemical Society

through the shape of the adsorbate alone. Subsequent studies on anthracene-based adsorbates have explored their use as molecular carriers,9 tested whether chemical modification of one end of the molecule can bias the ratio of forward to backward steps,10 and compared the diffusion constants of molecules with two and four thiol groups.11 The multitude of possible combinations of molecular shape and choice of surface structure point to the need for a more efficient computational strategy for searching the relevant parameter space than that provided by ab initio methods. Previously, Fichthorn and co-workers have used classical simulations to study the dynamics of linear alkanes physisorbed on metal surfaces.12,13 They found that concerted atom motions resulted in a lowering of the barrier heights for diffusion and noted an anisotropy in the diffusion constant. The flexibility of the alkanes tends to relax the constraints on dynamics imposed by a molecule’s shape and its registry with a crystal surface. In this paper, we set out to explore how rigid molecular shapes can be translated into directional “gaits” on a crystalline surface. We introduce a simple molecular model consisting of surface-binding groups or “feet” constrained by the rigid shape of the molecule. We consider a distribution of shapes, including a simple model of the DTA molecule. Our goal is to establish whether such a minimal model is sufficient to reproduce observed trends in mobility and, hence, provide an effective tool for understanding Received: November 14, 2010 Revised: March 29, 2011 Published: April 21, 2011 9526

dx.doi.org/10.1021/jp1108619 | J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C

ARTICLE

Figure 1. 9,10-Dithioacetylanthracene (DTA). Note that the acetyl groups (in gray) are removed after deposition on the Cu(111) surface in ref 4, and so it is the sulfurs that interact with the surface.

how adsorbed molecules move and how their shape determines the nature of this movement.

Figure 2. Model diamond-shaped adsorbate molecule (“walker”) with four interaction sites (“feet”) labeled 1
4, with positionsBri = (xi,yi) on the 2D surface. The sites are connected by rigid bonds and their positions specified by the coordinatesBr0 = (x0,y0) of the point labeled 0 and the angle θ between the y-axis and the vectorBr1
Br0. The length of the bond joining foot i to point 0 is li.

Table 1. Dimensions of Model Molecules

2. THE MODEL We model the adsorbates (“walkers”) as rigid structures decorated by sites (“feet”) that serve as the sole interaction with the surface. These feet are the points of contact between the molecule and the surface, and the manner of their collective movement over the surface determines the mobility of the adsorbed molecule. A foot might, for example, represent the sulfur atom via which an alkanethiol molecule binds to a metal surface. A pair of feet linked by a rigid spacer might also be used to model the relative rates of translational and rotational diffusion of a rigid feature such as the anthracene moiety in DTA. We neglect the flexibility of the molecule in the spirit of keeping the model as simple as possible. We note that flexibility in either the adsorbate or the substrate would have the general consequence of reducing the effect of anisotropies. The total potential energy of a walker on the surface is taken as V ðf B r i gÞ ¼ A0

4 X 3 X i

j

ci sinð B kj 3 B r iÞ

ð1Þ

where A0 sets the interaction energy scale while {ci} are dimensionless factors that determine the relative strength of the interaction of different feet with the surface. For the (111) surface of a face-centered-cubic (fcc) lattice, the unit √ vectors {k Bi} k 2 = (cos(π/3), sin(π /3)) are B k 1 = (1,0), B √ = (1/2, 3/2), and k 3 = (cos(
π/3), sin(
π /3)) = (1/2,
3/2). (Note that only B two of the three unit vectors are linearly independent, but the definition of three unit vectors on this two-dimensional surface allows the simple form of the potential energy in eq 1 to be used.) Working with these reduced vectors means scaling all lengths by the repeat length a along any one of these unit vectors. This repeat length (the 2D lattice spacing) is equal to the atomic spacing of the (111) surface, which is 2.55 Å in the case of copper. The rest of the parameters in the model can also be related to measurable physical quantities. Although eq 1 does not give the absolute surface binding energy, it does describe the positionor shape-dependent contributions to the binding energy of a molecule on the surface. In particular, the barrier to translational motion between adjacent minimum energy sites on the surface for a single “foot” with ci = 1 can be calculated and is 2.6A0 in the model. This energy barrier has been measured in density functional calculations of CH3S (i.e., a single “foot”) on a Cu(111) surface and was found to vary only slightly depending on whether

molecule I II III

l1 (= l3) √ 7/2 √ 7/2 √

l2 and l4 √ l2 = l4 = 7/2 l2 = l4 = 2

square diamond (commensurate)

l2 = 2, l4 = 5/2

diamond (incommensurate)

shape

(similar to DTA) 7/2

the C
S bond was constrained to be perpendicular to the surface (130 meV)14 or not constrained to any particular orientation (100 meV),15 although the absolute binding energy of CH3S does appear to be sensitive to the angle between the C
S bond and the metal surface.16 Using the barrier energy from the unconstrained calculation, the interaction amplitude in eq 1 can be estimated as A0 ≈ 40 meV (roughly 0.9 kcal/mol). A temperature of kBT = A0 therefore corresponds to about 490 K. All distances, energies, and times in this paper will be measured in units of the (111) repeat length a, A0, and τ0 = (m0a2/A0)1/2, where m0 is the effective mass of the adsorbate. For DTA on Cu(111), τ0 is estimated to be 1 ps, using a = 2.55 Å, taking A0 = 40 meV, and making the crude approximation that m0 is 1/4 of the molecular weight of DTA (distributing the weight equally among the 4 feet), which gives m0 = 60 amu. We shall consider three molecular shapes, labeled I, II, and III, each with four “feet” (tetrapods, in other words) as depicted in Figure 2. The three shapes differ in the length of a pair of “legs”. (The ith “leg” refers to the rigid spacer joining the ith “foot” to the center of the molecule.) The length of the legs for each molecule is provided in Table 1. The lengths of the first and third legs were chosen to mimic approximately the positioning of the thiol groups of DTA on the Cu(111) surface. As the overall role of flexibility in either the molecule or the surface is to diminish shape selectivity, we have left the consideration of such fluctuations to future papers. Model I represents the symmetric square tetrapod with legs of a single length. This length is chosen such that the pair of feet on one diagonal can simultaneously sit in local energy minima. Model II is the smallest bilateral (diamond) tetrapod with leg lengths chosen such that the pairs of feet at the end of one diagonal can simultaneously sit in minima. This model is distinguished from the square molecule of model I in that the two diagonals are of different length. In contrast, model III is an example of a bilateral molecule whose long axis (i.e., the axes along legs 2 and 4) is incommensurate with the structure of the solid surface; i.e., feet 2 and 4 cannot both sit in energy minima at the same time. 9527

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C

ARTICLE

Figure 3. (a, b) Two degenerate global minimum energy configurations (each 6-fold degenerate) for the type II1.5 walker with θ = 0.558π and 0.442π radians, respectively. (c) Global minimum-energy configuration (6-fold degenerate) for the type II2.5 walker with θ = 0.5π radians. Contours depicting the potential energy of interaction with the surface of a single foot with c1 = 1.0 are also shown (contour interval A0).

To complete the model, we must specify the choices of the individual “foot” interaction parameters {ci}. If all four “feet” are thiol groups, then a reasonable choice would be to set all the ci’s equal to one. We are interested, however, in using this model to explore a broader range of molecules such as, for example, DTA. As shown in Figure 1, the DTA molecule has only two thiol “feet”. It would be interesting, however, to see if the anthracene group can be usefully represented as a pair of feet with a different interaction strength with the surface. To this end, we fix c1 = c3 = 1.0 and explore the behavior of the walkers for the following values for c2 = c4: 0.5, 1.0, 1.5, 2.0, and 2.5 for each of the three types of walker. The value of c2 = c4 for a particular walker will be indicated as a subscript following its type (e.g., a I0.5 walker is a type I walker with c2 = c4 = 0.5). For each walker, all local √ minima in the √ unit cell of the√surface bound by √ the points (1/ 3, 1/2), (1/ 3, 3/2), (5/(2 3), 1), and (5/(2 3), 2) were found. A combination of steepest descent and conjugate gradient (Polak
Ribiere)17 minimization methods was employed, using analytical first derivatives of the potential in eq 1 with respect to x0, y0, and θ. The potential energy surface was searched exhaustively for local minima by running a series of minimizations with starting positions chosen randomly from a uniform distribution in the unit cell. For all of the walkers studied except for the II2.0 and II2.5 walkers, the global minimum-energy configurations have feet 1

and 3 approximately occupying potential energy minima on the lattice and have the long axis of the walker oriented at a slight angle to the high-symmetry axes of the surface, such as in the configurations depicted in Figure 3a,b for the II1.5 walker. These global minima come in degenerate or almost degenerate pairs (each of which is 6-fold degenerate, due to the 3-fold rotational symmetry of the surface and the reflection symmetry of the walkers) that can be approximately interconverted by a small angular rotation of around 0.1π radians about either foot 1 or foot 3 (i.e., a bipedal “walking” motion of feet 1 and 3 reminiscent of the way in which Kwon et al.4 postulated DTA moves on Cu(111)). For the II2.0 and II2.5 walkers, configurations like those shown in Figure 3a,b are local minima rather than global minima. For the global minima, feet 2 and 4 occupy potential energy minima on the surface and the long axis of the walker is aligned with the highsymmetry axis of the lattice, while feet 1 and 3 no longer occupy potential energy minima on the lattice, as shown in Figure 3c. For the other walkers, configurations like that shown in Figure 3c are local minima. These results match results obtained by Kwon et al.4 for DTA on Cu(111). Their experiments suggest a global minimum-energy configuration of DTA on Cu(111) like that depicted in Figure 3c while DFT calculations that they carried out for the DTA/Cu(111) system give global minima like those depicted in Figure 3a,b. 9528

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C

ARTICLE

A detailed study of the dependence of the ground-state arrangement of the molecule and its energy will be presented elsewhere (work in preparation). As a simple and computationally efficient means of simulating the motion of the walkers on the surface, Brownian dynamics with a delta-correlated stochastic force was employed. The feet of the walkers were moved independently according to the Langevin equation18 mi v_ i ðtÞ ¼
mi γi vi ðtÞ þ fi ð½qi ðtÞÞ þ Ri ðtÞ

diffusion tensor, D, for the center-of-mass of the walker, 4 4 ir i(t)/∑i mi Br c = ∑i mB " # Dxx Dxy D¼ ð3Þ Dyx Dyy where ÆΔζðtÞΔωðtÞæ ¼ 2Dζω t,

ð2Þ

where the index i = 1, 2, ..., 8 labels the x and y coordinates of the feet, and mi, vi = q_ i, and γi are the mass, velocity, and friction coefficients corresponding to the coordinate qi. The deterministic force fi(q B) is the gradient of the potential in eq 1 with respect to the coordinates and the stochastic force Ri(t) has zero mean and a delta-function correlation in time, ÆRi(t)Rj(t0 )æ = 2mikBTγiδijδ(t
t0 ), where kB is Boltzmann’s constant and T is the temperature. Constraints on the bond lengths and angles of the rigid walkers were enforced by using the SHAKE algorithm.19 The third-order Brownian dynamics algorithm due to van Gunsteren and Berendsen,20 which incorporates SHAKE, was used to integrate the equations of motion. For all the simulations, the masses of the feet and friction coefficients on the feet were taken to be equal to each other, i.e., mi = m0 and γi = γ = 72.6/τ0. For the potential energy surface used (eq 1), this choice of γ is consistent with experimentally fitted friction coefficients used in previous Langevin dynamics studies of surface diffusion of small molecules21,22 and puts the system in a regime of relatively high friction (strong coupling to the thermal bath). The choice of γ can potentially affect the mechanisms and activation barriers for diffusion. For example, it can influence whether the molecules exhibit single or multiple unit-cell hops on the surface, which in turn impacts the Arrhenius parameters for activated diffusion.23 However, the large potential energy barriers to diffusion of the chemisorbed molecules that we focus on in this article (>4kBT at 298 K for alkanethiol “feet” on Cu(111) surfaces) make multiple hops unlikely. A complete characterization of the effect of molecular shape on surface diffusion (for example, to treat more weakly physisorbed molecules or high temperatures) would account for the low-friction regime and the possibility of multiple hops. In the high-friction regime, surface diffusion coefficients should simply scale with 1/γ. In all the simulations, the walker’s feet were given identical initial velocities (i.e., the walker initially had zero angular momentum), which were chosen from a Maxwell
Boltzmann distribution for the specified temperature T. This arbitrary choice of initial velocities did not affect the dynamics significantly because the system was in the high-friction regime, and so the relaxation of the velocity occurred on a much faster time scale than any relevant motion of the walker across the surface. For each walker, trajectories were started from three of the global minimum-energy configurations that could be interconverted by a 2π/3 radian rotation of the surface (i.e., the three different starting configurations had their long axes aligned (approximately) along each of the three high-symmetry axes of the lattice). For each starting configuration, a set of 1000 independent trajectories were initiated, from which correlation functions were calculated. For each of the starting configurations, we computed the elements of the translational

ζ, ω ¼ x or y of the center of mass 0

ð4Þ

0

and Δζ(t) = ζ(t þt)
ζ(t ). The duration of the simulations was chosen such that the mean-squared displacement of the center of mass of the walker at the end of the simulation was approximately 2a2 (i.e., the walker moved on average more than a unit cell over the length of the simulation). The diffusion coefficients, Dþ and D
, along the pair of orthogonal directions for which the rate of translational diffusion was maximum and minimum, respectively, were calculated as the eigenvalues of the diffusion tensor D in eq 3 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ð5Þ D( ¼ ½Dxx þ Dyy ( 4Dxy 2 þ ðDxx
Dyy Þ2  2 The total diffusion coefficient, D = (Dxxþ Dyy)/2, was also calculated. The rotational relaxation time τr was calculated by 0 u (t0 þt)æ fitting the angular autocorrelation function Æu B(t ) 3 B (where the unit vector B u points from the molecule center to foot 1) to the function exp(
t/τr). For each walker, the aforementioned simulations were carried out at 4 or 5 different temperatures between kBT = 0.6A0 and 4.0A0, with the particular range chosen for each walker such that activated (as opposed to free) diffusion was observed at all temperature and such that the temperature points were roughly evenly spaced in 1/T.

3. RESULTS 3.1. Anisotropic Motion: Shape vs Interactions. Within our simple reduction of the molecule
surface interaction to a discrete number of feet, the anisotropy of this interaction can be classified two ways. The shape anisotropy of the molecule refers to the physical distribution of the feet around the center of mass. For our model tetrapods we shall measure this by the ratio of the length of the long axis over that of the short. The interaction anisotropy refers to the difference between the interactions of the various feet and the surface. Given the constraint we have imposed on the interaction strength c1 between foot 1 and the surface (see section 2), the interaction anisotropy can be measured as the deviation of c2 from one. Over the range of temperatures studied, the translational diffusion constant D, rotational diffusion constant R = 1/τr, and diffusion anisotropy Dþ/D
are well described by an Arrhenius expression of the form ν exp(
Ea/kBT). The values of the frequency factor ν and the activation energy Ea for translational and rotational motion for the three walkers with different values of c2 are provided in Figures 4a and 4b, respectively. We find that the activation energy for translational and rotational motion are the same for molecules I and III regardless of the interaction anisotropy, while rotational motion has a higher activation energy than translation for molecule II, with the difference increasing with interaction anisotropy. The frequency factor is smaller for rotation than translation in all cases. This means that, even with the same activation energy, molecules I 9529

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C

ARTICLE

Figure 6. Activation energy Ea for translational diffusion vs the global minimum energy Emin of the walker on the surface for the type I, II, and III walkers.

Figure 4. (a) Natural logarithm of frequency factor, ln(ν), and (b) activation energy, Ea, for translational and rotational diffusion vs the interaction strength, c2 = c4, of feet 2 and 4 with the surface for the type I, II, and III walkers. The error bars are one standard deviation of the errors in the coefficients from linear regression of Arrhenius plots for translation and rotation, respectively.

Figure 5. Difference between the activation energies, Δ = Ea

Eaþ, for translation along the orthogonal directions of maximum and minimum diffusion vs the interaction strength, c2 = c4, of feet 2 and 4 with the surface for the type I, II, and III walkers. The error bars are one standard deviation of the errors in the coefficients from linear regression of Arrhenius plots of D
/Dþ.

and III can still exhibit a rotational relaxation that is significantly slower than translational motion and, hence, anisotropic translational diffusion. The anisotropy of the diffusive motion can be quantified by the temperature dependence of the ratio D
/Dþ, i.e., the difference between the activation energies Δ = Ea

Eaþ and the ratio of the frequency factors. In Figure 5 we plot Δ against c2 for the three walkers. Varying c2 from 1.0 corresponds to increasing the interaction anisotropy, and we find an increase in the difference Δ between the activation energies between the slow and fast directions with increasing c2. Walker II exhibits the largest increase and hence the greatest anisotropy in its motion. The ratio of frequency factors ν
/νþ (not shown) exhibits a similar dependence with c2. The mechanism by which this selection of direction is achieved will be examined below.

One aspect of the data in Figure 5 requires some consideration. Increasing c2 not only increases the difference between the activation energies along the slow and fast directions but also increases the magnitude of the activation energies overall. Increasing directionality is of limited interest if the mobilities one refers to are all steadily going to zero. In Figure 4b, we plot the activation energy of the (total) translational diffusion constant D against c2. While the walkers I and III exhibit the expected increase in activation energy with increasing interaction of feet 2 and 4 with the surface, walker II, the commensurate diamond, exhibits an activation energy that is virtually independent of c2. This independence is remarkable. It means that walker II, in stark contrast to expectations and the behavior of the other molecule models studied here, has an activation energy for translational diffusion that is only very weakly dependent on its adsorption energy (defined here as
Emin where Emin is the global energy minimum), as shown in Figure 6. The increasing anisotropy of the interactions of the feet in walker II increases the anisotropy of the diffusive motion without changing the overall mobility; the decrease in D
is compensated by the roughly equivalent increase in Dþ. Perhaps the clearest measure of this behavior is to look at the ratio Dþ/D
as a function of c2 where the temperature is adjusted for each point so that the total diffusion constant is constrained to have the same value for all walkers (7.5  10
5a2/τ0). Shown in Figure 7, this plot reveals the dramatic increase in the anisotropy of the motion of walker II, in stark contrast with the absence of any significant increase in the diffusion anisotropy of walkers I and III with increasing c2. Our results highlight how the two aspects of molecular anisotropy—shape and interaction—influence an adsorbate’s mobility. Given the right shape, an increase in the interaction anisotropy produces a corresponding increase in anisotropic mobility. Walker II meets this requirement of the “right” shape while walkers I and III do not. Without interaction anisotropy (i.e., c2 = 1.0), walker II manages only a modest anisotropy Dþ/D
∼ 2. We conclude that neither shape or interaction anisotropy can, on their own, produce the significant anisotropy that they can achieve when properly combined. The results presented so far were for simulations in which the molecules moved on average slightly more than a unit cell over the length of the simulation. The measured diffusion coefficients therefore characterize the anisotropy of an elementary hop between adjacent unit cells. The diffusion anisotropy is expected 9530

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C

Figure 7. Diffusion anisotropy, Dþ/D
, at the temperature at which the total diffusion constant is D = 7.5  10
5a2/τ0 vs the interaction strength, c2 = c4, of feet 2 and 4 with the surface for the type I, II, and III walkers. The points were calculated using the fitted Arrhenius coefficients.

to diminish with subsequent hops, as the angular orientation of the diffusing molecule becomes uncorrelated with its initial orientation, and after long enough times the diffusion anisotropy will tend to 1. But the large diffusion anisotropies and small rotational diffusion coefficients for a single hop of some of the molecules studied indicate that this decorrelation can be slow. We have rerun the simulations of the II1.5, II2.0, and II2.5 walkers (the molecules that exhibited the highest diffusion anisotropies), so that the average mean-squared displacement over the simulations was 9a2, i.e., the molecule moved on average 3 unit cells. For the conditions shown in Figure 7, the measured diffusion anisotropies were 4.9, 7.8, and 8.5 respectively for the II1.5, II2.0, and II2.5 walkers in the longer simulations. As expected, these values are lower than the values of 6.4, 15.4, and 15.1, respectively, shown in Figure 7 for the shorter simulations. But, even after moving 3 unit cells, the diffusion anisotropies of these molecules are still substantial. We have shown that our simple model for the type II walkers qualitatively reproduces the experimentally observed diffusion anisotropy of DTA on a Cu(111) surface.4 The model also demonstrates reasonable quantitative accuracy. If the activation energies and frequency factors for translation are converted from reduced units, in which they are roughly 4A0 and e
3.7a2/τ0, respectively, for all type II walkers, to real units using the conversion factors given in section 2, 160 meV and 2.5  1010a2/s are obtained, respectively. Given the simplicity of the model, these values compare well with the measured values of 130 meV and 4  109a2/s for DTA.4 The agreement between the simulations and experimental results indicates that the chosen simulation parameters describe the real physical system well; in particular, the agreement for the frequency factor justifies our assumption of dynamics in the high-friction regime. 3.2. Effect of an Applied Field. Using our model, we can also predict how the motion of molecules on surfaces is modified by an applied field and how the field couples to the molecular anisotropy. This applied field could be an electric field or chemical potential gradient arising from inhomogeneous adsorbate coverage of the surface.24 In particular, physically reasonable applied electric fields have been shown in ab initio calculations to significantly affect the energy barriers to diffusion of alkanethiols (i.e., molecules comprising single “feet” in our model) on metal

ARTICLE

Figure 8. Diffusion constant in the x-direction, Dxx, and in the y-direction, Dyy, in the presence of an external field in the x-direction, φx, or in the y-direction, φy, for type (a) II2.5 and (b) III2.5 walkers. Trajectories were started from the global minimum-energy configuration with θ = 0.5π radians for the II2.5 walker and θ = 1.105π radians for the III2.5 walker. (Diffusion constants for the strongest fields are approximate due to deviations of the mean-squared displacement from linearity versus time; otherwise, errors in the diffusion constants were estimated to be less than 10%.)

surfaces.16 We examined the effect on the dynamics of the II2.5 and III2.5 walkers of adding a linear field, which adds a term ∑4i (φxxi þ φyyi) to the potential energy in eq 1. Note that this field acts equally on all feet. Elements of the diffusion tensor were obtained as before, except that simulations were started only from the global minimum energy configurations with θ = 0.5π radians and 1.105π radians respectively for the II2.5 and III2.5 walkers. These starting configurations led to translational diffusion predominantly in the y-direction in the field-free simulations. Fields were applied in either only the x-direction or only the y-direction. Figure 8 shows the diffusion coefficient in the x-direction, Dxx, and in the y-direction, Dyy, as a function of the field strength for the two walkers. Not surprisingly, applying a field in the x-direction results in an increase in diffusion in the x-direction for both walkers. Similarly, applying a field in the y-direction results in an increase in diffusion in the y-direction. More surprisingly, for the II2.5 walker, applying a field that increases in the positive x-direction results in large increases in diffusion in both the x- and y-directions, with the increase in diffusion in the y-direction the larger of the two. On the other hand, applying a field that increases in the negative x-direction results in a decrease in diffusion in the y-direction at the same time as an increase in diffusion in the (negative) xdirection. For the III2.5 walker, applying a field in the x-direction has little effect on diffusion in the y-direction. For both walkers, applying a field in the y-direction has little effect on diffusion in the x-direction. The dependence of the diffusion of the II2.5 walker in the y-direction on the applied field in the x-direction is reminiscent of the current
voltage curve of a diode, with the diffusion rate exhibiting a highly asymmetric dependence on the field direction. This unusual behavior could potentially be exploited as a switch for turning on and off molecular motion. Given the apparent 9531

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C

ARTICLE

Figure 9. Typical trajectory of (x0,y0) starting from the global minimum energy configuration with (a) θ = 0.5π radians for the II2.5 walker and (b) θ = 1.105π radians for the III2.5 walker. Contours depicting the potential energy of interaction with the surface of a single foot with c1= 1.0 are also shown (contour interval A0).

similarity of the motion of the II2.5 walker and DTA on a (111) surface, it would be interesting to see whether the surface diffusion of DTA molecules can be controlled using this effect. Other unusual geometry-dependent effects of applied fields on adsorbate diffusion may be possible. 3.3. Molecular Gaits: Coupling Rotations to Translations. To explain the peculiar effects of an applied field on the diffusion of the II2.5 walker, as well as the diffusion anisotropy observed in field-free simulations, we have examine the detailed mechanisms for adsorbate motion. Figures 9a and 9b show typical trajectories for the II2.5 and III2.5 walkers, respectively. The initial orientation of the walkers in the two figures was θ = 0.5π radians and 1.105 π radians, respectively. At first glance, the trajectories for the two walkers seem quite similar, although the III2.5 walker shows less unidirectionality in its motion. However, the initial orientation (which is to a large degree maintained throughout the trajectory) differs by around π/2 radians for the two walkers: the II2.5 walker moves most rapidly in the direction of its long axis—one of the [110] directions (much like DTA on Cu(111))—while the III2.5 walker moves most rapidly in the direction of its short axis. The mechanism of motion therefore appears to be quite different for the two walkers. Another difference is that the global minimum for the II2.5 walker is positioned to one side of the path along which the molecule diffuses. This walker spends most of its time in these lateral “dead ends” (its global minimum) and has to move in a direction orthogonal to the general direction of translational diffusion to regain the path along which it can move. In contrast, the global minimum for the III2.5 lies along the path of translational diffusion. In order to understand better the differences in the diffusive behavior of the II2.5 and III2.5 walkers, we computed all the minimum-energy pathways between minima connected by single saddle points on the potential energy surface for the walkers. To do this, we determined the positions in (x0,y0,θ)-space of all of the saddle points in a unit cell of the lattice by minimizing the square of the absolute value of the gradient of the potential energy in eq 1, |rV|2. (This procedure yields all the minima and maxima as well as the saddle points.) From each saddle point (ignoring symmetry-related saddle points, which are easily related by rotation or reflection), we ran Brownian dynamics trajectories at close to zero temperature and moderate friction to obtain the minimum-energy pathways.

Figure 10. Mechanism for diffusion between adjacent unit cells that requires surmounting the lowest energy barriers for (a) II2.5 and (b) III2.5 walkers starting from global minimum-energy configurations with θ = 0.5π radians and 1.105π radians, respectively. The global minimum energy configurations in adjacent unit cells are drawn as solid circles connected by solid lines, while the configuration corresponding to the highest barrier along the diffusion path is drawn as empty circles connected by dashed lines. The crosses are the positions of the potential energy minima for a single “foot”.

Figure 10 illustrates the mechanisms for diffusion of the II2.5 and III2.5 walkers between adjacent unit cells on the surface that require surmounting barriers of the lowest energy, starting from the same global minimum-energy configurations that were the starting configurations of the trajectories in Figure 9. In both cases, these pathways involve only small deviations in the orientation θ and connect global minima in adjacent unit cells of the same or only slightly different orientation in the case of the III2.5 walker. Maintaining the same orientation with respect to the surface is an essential component of anisotropic diffusion. Large 9532

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C

Figure 11. Lowest energy pathway to translation shown in Figure 10 (solid line) and minimum-energy pathway that converts the global minimum energy configuration to a degenerate global minimum that differs in orientation by a rotation of 2π/3 radians (dashed line) for (a) II2.5 and (b) III2.5 walkers. Minima and barriers for the lowest energy pathways are labeled with “M” and “B”, respectively. The numbers following the M or B labels correspond to the order on the potential energy surface (from lowest to highest energy) of the minimum or barrier (e.g., M1 is the global minimum and B1 is the lowest energy barrier). B3 in (a) and B6 in (b) are depicted in Figures 10a and 10b, respectively.

rotations allow the walker to select another symmetry-related direction for translational motion. In both cases, the pathway involves motion generally along the y-direction, as expected from Figure 9. But, as pointed out above, the II2.5 walker moves in the direction of its long axis, whereas the III2.5 walker moves in the direction of its short axis. Kwon et al.4 postulated that DTA moves forward along the surface by what is essentially a bipedal walking gait in which the thiol moieties act like feet. Unfavorable energetic interactions of the anthracene moiety with the surface prevent the molecule from rotating too much between each step and therefore the molecule moves exclusively in the direction of the anthracene alignment. Molecular gaits resulting from thermal fluctuations result in stochastic locomotion. The observation of regular periodic behavior in biological walkers like kinesin25 is a consequence of coupling motion with an effectively irreversible dephosphorylation. The gaits of both II2.5 and III2.5, as shown in Figure 10, consist of a translation coupled to small periodic oscillations of the molecular orientation, i.e., a bipedal walk. For the II2.5 walker, a “dead end detour” (associated with the global minimum in Figure 10) is reached via a coupled rotation and lateral displacement. Escape from this minimum is rate determining with regards to the mobility along the y direction. The III2.5 walker suffers a similar transient trapping, the difference being that for the III2.5 walker there is little lateral motion associated with the trap. Figure 11 depicts the energy changes that occur along the pathways in Figure 10 (which for II2.5 and III2.5 join global minimum-energy configurations in adjacent unit cells that have the same orientation) as well as minimum-energy pathways that connect global minimum energy configurations with different

ARTICLE

orientations (not depicted in Figure 10). The barriers along the latter pathway measure how readily the initial orientation of the walker is randomized. The difference between the barriers along the two pathways therefore measures the degree of anisotropy of diffusion. Figure 11 shows that the lowest energy pathway between adjacent unit cells for both walkers involves translation without rotation. But while this pathway requires surmounting only the three lowest energy barriers on the potential energy surface for the II2.5 walker (i.e., motion along any other path on the potential energy surface involves higher barriers), it involves the second, fourth, and sixth lowest barriers for the III2.5 walker. Figure 11 also shows that the barrier that the walker must surmount to change orientation is significantly higher than that required to undergo translational diffusion without rotation for the II2.5 walker; the difference in barrier heights is substantially smaller for the III2.5 walker. The difference in energy between the highest barriers along the two pathways for the two walkers are also consistent with the activation energy difference Δ between the slow and fast diffusion directions depicted in Figure 5. Furthermore, the highest barriers on the lowest energy pathways between adjacent unit cells for the II2.5 and III2.5 walkers are within 15% of the activation energies for translational diffusion in Figure 4b, indicating that diffusion occurs by hopping between binding sites in which the walkers are mostly localized, as suggested by the trajectories in Figure 9. Larger deviations would be expected for nonlocalized diffusion.12,23 An examination of the minimum energy pathways in Figure 10 or the trajectories in Figure 9 explains why diffusion of the II2.5 walker in the y-direction is enhanced by applying a positive field in the x-direction: the field pushes the walker out of the global minimum in which it spends most of its time and toward the path in which it moves to translate in the y-direction.

4. CONCLUSION In this paper we have shown how a molecule’s shape can dramatically influence both the magnitude and anisotropy of its mobility on a high-symmetry surface. The interaction between the structures of the adsorbed “walker” and the surface can result in a complex range of behaviors as either the shape or the interaction anisotropy is varied. We find that the shape dominates the observed behavior. Anisotropy in the adsorbate
surface interaction strengths can significantly augment the effects of shape anisotropy but they cannot replace it. The manifestation of persistent translation along a unique direction on a (111) surface requires the fluctuations in the molecular orientation to be restricted to a range smaller than π/3. We have shown that a diamond-like shape, similar to that of DTA, accomplishes this. The resulting mobility of this particular shape (our II2.5 walker) lies along a linear path corresponding to the [110] direction on the (111) surface. We find that a field applied perpendicular to the preferred direction of translational diffusion can accelerate or retard the mobility of the molecule along the preferred direction in cases where the rate-determining step involves short perpendicular “detours” from this path. In addition to the appearance of directional selection through the interaction of the molecular shape with the crystal surface, we have also demonstrated that suitably shaped molecules can accomplish the remarkable feat of (approximately) decoupling the activation energy related to the motion from their binding energy to the surface. The assumption that a strongly bound molecule must, unavoidably, be laterally pinned is, therefore, 9533

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534

The Journal of Physical Chemistry C incorrect. A molecule, characterized by multiple binding sites (referred to as “feet” in this paper), can be given a shape such that some feet are sliding into one energy minimum as other feet are climbing out of another, thus offsetting the energy cost of the transition and facilitating mobility. We have also found that our model reproduces both qualitatively and semiquantitatively the diffusive motion of the molecule DTA on a Cu(111) surface. The simplicity of the model means that all possible pathways for diffusion on the surface can be enumerated with ease, allowing a precise determination of the mechanism by which a molecule “walks” on the surface and how this mechanism depends on molecular shape. A detailed understanding of the mechanism for adsorbate diffusion makes it possible to design methods to control diffusion, for example, by applying an external field along a particular axis, as illustrated in this work.

ARTICLE

(17) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, P. B. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, 1992; pp 397—408, 420—425. (18) Allen; M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon: Oxford, 1987; pp 259—264. (19) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327–341. (20) van Gunsteren, W. F.; Berendsen, H. J. C. Mol. Phys. 1982, 45, 637–647. (21) Ferrando, R.; Spadacini, R.; Tommei, G. E. Surf. Sci. 1992, 269/ 270, 184–188. (22) Graham, A. P.; Hofmann, F.; Toennies, J. P.; Chen, L. Y.; Ying, S. C. Phys. Rev. Lett. 1997, 78, 3900–3903. (23) Raut, J. S.; Fichthorn, K. A. J. Chem. Phys. 1995, 103 8694–8704. (24) Tsong, T. T. Physica A 2005, 357, 250–281. (25) Lipowsky, R.; Beeg, J.; Dimova, R.; Klumpp, S.; M€uller, M. J. I. Physica E 2010, 42, 649–661.

’ AUTHOR INFORMATION Corresponding Author

*Phone: þ61-8-8303-5580. Fax: þ61-8-8303-4380. E-mail: [email protected].

’ ACKNOWLEDGMENT We thank David Barrett, Nicholas Lambropoulos, and Jeff Reimers for helpful discussions. ’ REFERENCES (1) Vericat, C.; Vela, M. E.; Benitez, G.; Carro, P.; Salvarezza, R. C. Chem. Soc. Rev. 2010, 39, 1805–1834. (2) Love, J. C.; Estroff, L. A.; Kriebel, J. K.; Nuzzo, R. G.; Whitesides, G. M. Chem. Rev. 2005, 105, 1103–1169. (3) Lambropoulos, N. A.; Reimers, J. R.; Hush, N. S. Nanotechnology 2004, 15, 1226–1232. (4) Kwon, K.-Y.; Wong, K. L.; Pawin, G.; Bartels, L.; Stolbov, S.; Rahman, T. S. Phys. Rev. Lett. 2005, 95, 166101. (5) Romero, A. H.; Lacasta, A. M.; Sancho, J. M. Phys. Rev. E 2004, 69, 051105. (6) Otero, R.; H€ummelink, F.; Sato, F.; Legoas, S. B.; Thostrup, P.; Læsgaard, E.; Stensgaard, I.; Galv~ao, D. S.; Besenbacher, F. Nature Mater. 2004, 3, 779–782. (7) Schunack, M.; Rosei, F.; Naitoh, Y.; Jiang, P.; Gourdon, A.; Læsgaard, E.; Stensgaard, I.; Joachim, C.; Besenbacher, F. J. Chem. Phys. 2002, 117, 6259–6265. (8) Stranick, S. S.; Kamna, M. M.; Weiss, P. S. Science 1994, 266, 99–102. (9) Wong, K. L.; Pawin., G.; Kwon, K.-Y.; Lin, X.; Jiao, T.; Solanki, U.; Fawcett, R. H. J.; Bartels, L. Science 2007, 315, 1391–1393. (10) Pawin, G.; Wong, K. L.; Kwon, K.-Y.; Frisbee, R. J.; Rahman, T. S.; Bartels, L. J. Am. Chem. Soc. 2008, 130, 15244–15245. (11) Cheng, Z.; Chu, E. S.; Sun, D.; Kim, D.; Zhu, Y.; Luo, M.; Pawin, G.; Wong, K. L.; Kwon, K.-Y.; Carp, R.; Marsella, M.; Bartels, L. J. Am. Chem. Soc. 2010, 132, 13578–13581. (12) Huang, D.; Chen., Y.; Fichthorn, K. A. J. Chem. Phys. 1994, 101, 11021–11030. (13) Raut, J. S.; Fichthorn, K. A. J. Chem. Phys. 1998, 108 1626–1635. (14) Ferral, A.; Patrito, E. M.; Paredes-Olivera, P. J. Phys. Chem. B 2006, 110, 17050–17062. (15) Kon^opka, M.; Rousseau, R.; Stich, I.; Marx, D. Phys. Rev. Lett. 2005, 95, 096102. (16) Cometto, F. P.; Paredes-Olivera, P.; Macagno, V. A.; Patrito, E. M. J. Phys. Chem. B 2005, 109, 21737–21748. 9534

dx.doi.org/10.1021/jp1108619 |J. Phys. Chem. C 2011, 115, 9526–9534