Modeling the Diffusive Motion of Large Organic Molecules on

Nov 16, 2008 - Q. Guo , A. Paulheim , and M. Sokolowski , H. Aldahhak , E. Rauls , and ... Florian Schiffmann , Joost VandeVondele , Jürg Hutter , Ro...
0 downloads 0 Views 1MB Size
Download PDF
J. Phys. Chem. C 2008, 112, 19577–19583

19577

Modeling the Diffusive Motion of Large Organic Molecules on Insulating Surfaces Thomas Trevethan* and Alexander L. Shluger London Centre for Nanotechnology, UniVersity College London, Gower Street, London, U.K. ReceiVed: July 18, 2008; ReVised Manuscript ReceiVed: October 6, 2008

We show that the overall structure and flexibility of an organic molecule has a profound effect on the mechanism of diffusion and the effective diffusion rate on a surface. Calculations were performed to model the diffusion of a set of large organic molecules with polar binding groups on the perfect TiO2 (110) surface. These simulations involved determining the accessible states of the molecule surface system and the energy barriers that separate them using realistic atomistic simulations employing a set of specifically developed potentials. With the complete set of accessible states and energy barriers for each molecule, we then performed kinetic Monte Carlo simulations to determine the mechanisms of diffusion and the effective diffusion rates of the molecules. These calculations suggest ways in which the mobility of large molecules on surfaces can be controlled through careful design of the molecular structure. Introduction Understanding the adsorption and diffusion of large organic molecules on surfaces is important to many areas of surface science and nanotechnology. This is particularly true in the formation of self-assembled monolayers, where the rates of various diffusive processes may have a profound impact on the eventual structures created, and is also important in the areas of catalysis, coatings, and molecular electronics. The development of single molecule devices for molecular electronics1–3 will be reliant upon exact control over the position and mobility of large molecules on insulating surfaces, where it is often important that the molecule be bound and immobile for long periods of time, often at room temperature. In applications such as these, it is desirable to develop strategies for designing or choosing molecular structures to tailor their mobility on a particular surface. To do this, it will be necessary to accurately evaluate the diffusion of large and complex molecules on a specific surface with atomistic simulations, which is a very challenging task because of the complexity of this type of system and the timescales over which molecular motion occurs. Long time-scale simulation methods must be used: molecular dynamics simulations precisely follow the real-time trajectory of the system;4 however, they are only computationally feasible over picosecond to nanosecond time scales. Processes related to surface diffusion in the type of systems described above may occur over seconds, minutes, or hours. In this paper, we evaluate the diffusive motion of a class of large organic molecules on an oxide surface using a combination of atomistic modeling and kinetic Monte Carlo calculations. The diffusive motion of a molecule on a surface consists of transitions among many different minima on the complex potential energy surface of the surface-molecule system. For each of the molecules investigated here, these accessible states and the energy barriers that separate them are determined, and then the resultant finite temperature dynamical behavior of the molecule is evaluated with the kinetic Monte Carlo method. As a result of these simulations, we are able to show how the structure of the molecule is critical to its diffusive motion and how making small adjustments to the molecular structure can * E-mail: [email protected].

have a dramatic effect on the overall behavior of the system. We are able to suggest general ways in which molecular structures could be designed to either maximize or minimize the mobility of the molecule on a surface and also how molecular structures could be designed to move in only a single direction. To investigate the general issues relevant to this type of system, we have employed a model system that is formed of a set of closely related molecules, consisting of polycyclic aromatic “backbones” with four attached carboxylic acid based groups, adsorbed on the TiO2 rutile (110) surface. This substrate is chosen because it is technologically relevant and reasonably well understood5 and its structure reduces the diffusion of carboxylic acid groups on the surface to one dimension, which significantly simplifies the treatment of this system. The rutile surface structure is highly anisotropic: there are rows of Ti atoms and bridging O atoms along the direction. It is known that a carboxylic acid group will bind to the surface in a bidentate configuration, with the two O atoms of the acid group bonding to two adjacent surface Ti atoms along the Ti rows.6,7 At room temperature, diffusion of this group along the rows can occur, but the barrier for diffusion across the bridging oxygen rows is insurmountable,8 leading to 1-D diffusive motion. As is the case with most oxide surfaces, the interaction of nonpolar hydrocarbon molecules is very weak, and therefore, polar molecular groups are often necessary to bind a molecule to the surface.9,10 The molecules employed in this study are shown in Figure 1. Molecules A and B consist of an anthracene backbone with four carboxylic acid groups substituted at the 1, 2, 3, and 4 positions indicated. For molecule A, the substituted groups are -CH2COOH, and for molecule B, are -(CH2)2COOH. Molecule C consists of a tetracene backbone, with four -CH2COOH groups again substituted at the 1, 2, 3, and 4 positions of the molecule. These molecules were chosen to be representative of rigid polycyclic aromatic species that employ polar groups that act as “legs” to bind to a nonreactive oxide surface. Molecules A and B have the same backbone structure but differ in the length, and therefore the “flexibility”, of the legs. Molecules A and C have the same structure legs, but different length backbones, varying the distance between the pairs of legs.

10.1021/jp806355m CCC: $40.75  2008 American Chemical Society Published on Web 11/16/2008

19578 J. Phys. Chem. C, Vol. 112, No. 49, 2008

Trevethan and Shluger

Figure 1. Schematic representation of the three molecules.

Atomistic Calculations The atomistic calculations presented here have been carried out using a classical potential model that has been developed to describe the interaction of organic species with the perfect TiO2(110) surface. Intramolecular interactions are described using the AMBER force field,11 and the interactions between atoms in the surface are described using the force field of Bandura et al.12 Interactions among atoms in the molecule and atoms in the surface are described using a specially constructed force field for the interaction of carboxylic acid groups, saturated hydrocarbon chains, and phenyl groups with the rutile surface. This was developed by fitting pair potentials to potential energy surfaces obtained from ab initio calculations, which are described in detail in ref 13. The TiO2(110) surface was represented by a finite cluster nine atomic layers deep and providing a surface area of 4.78 nm × 5.24 nm (16 × 8 surface unit cells). This cluster was constructed such that it was symmetrical and charge-neutral.14 The upper five atomic layers of the central 8 × 4 unit cells were free to relax, and all other atoms in the cluster were frozen in bulklike positions. All of the calculations have been performed with the Sci-Fi atomistic simulation code,15 which relaxes the coordinates of all free atoms using the BFGS minimization algorithm. We first consider the adsorption of molecule A, where the lowest energy configuration is shown in Figure 2. Here, the two carboxylic binding groups at positions 1 and 2 on the molecule bind to two adjacent pairs of Ti atoms along a Ti row, and the other two binding groups at positions 3 and 4 bind to the

Figure 2. Ground state configuration of molecule A: (i) side view and (ii) end view .

corresponding two pairs of Ti atoms in the next Ti row, with the anthracene backbone straddling the oxygen row separating the two Ti rows, and parallel to it. An alternate stable adsorption configuration for this molecule (separated from the lowest energy configuration by a barrier that is insurmountable at room temperature) is with the 1 and 4 groups bonding to Ti atoms in

Diffusive Motion of Large Organic Molecules

J. Phys. Chem. C, Vol. 112, No. 49, 2008 19579

Figure 5. Minimum energy path for the transition of leg 1 from the ground state into the leading-leg state for both molecule A (solid line) and molecule B (dashed line).

Figure 3. Ground state configuration of molecule A in its rotated adsorption configuration: (i) side view and (ii) end view .

Figure 4. Top view of the ground-state configuration of molecule A. The Ti atoms in the surface directly bonded to the carboxylic acid binding groups are highlighted.

the same Ti row and the 2 and 3 groups in the adjacent row, with the anthracene backbone now perpendicular to the rows (Figure 3). This configuration is approximately 0.7 eV higher in energy than the ground state. Figure 4 shows an overhead view of molecule A in its lowest energy configuration, with the Ti atoms that bond directly to the binding groups highlighted. Since the hydrocarbon part of the molecule interacts very weakly with the surface and the carboxylic acid groups bind relatively strongly to the Ti surface atoms, the configuration of the molecule can be described in terms of the relative positions of the binding groups, which define a potential energy basin. As described in ref 8, a single carboxylic acid group can move along the Ti row to an adjacent state one lattice spacing away by crossing a moderate energy barrier. In the case of a single benzoic acid molecule, this barrier is approximately 0.6 eV. This results in the possibility that there exist different accessible potential energy basins (or states) for the molecule with the binding groups in different relative positions along the Ti rows. The existence of these states, their energies, and the energy barriers that separate them will depend on the structure and mechanical properties of the entire molecule. These states can be explored by following the minimum energy path for a

particular binding group to move into a neighboring position along the Ti rows, which can be achieved using a constrained minimization procedure, described in detail in ref 8. Figure 5 shows the minimum energy path for the transition of leg 1 of the molecule A to its neighboring position, moving left, away from the center of the molecule (this will be called the “leading leg” configuration). Because of the symmetry of the molecule, this potential energy surface is the same for each of the four legs to move in the same way. There is no accessible basin for any of the legs to move in the opposite direction in this configuration because of the presence of the other leg in the neighboring position along the Ti row. In addition, there is no accessible state for the leg to move a farther distance away from the center of the molecule. Once the molecule is in this new configuration, additional basins can be found by moving the other binding groups into new positions by following the minimum energy path. Using this procedure and the symmetry of the system, it is possible to map all of the accessible, nonequivalent states of the molecule surface system and the energy paths and, hence, barriers that separate each of them. At this point, we introduce a notation to simply describe the configuration of the molecule in terms of the relative positions of the binding groups. Each of the binding groups or legs (1, 2, 3, 4) is given an integer to describe its position, with 0 the initial position in the ground state, -1 to have moved one lattice distance left, and +1 to have moved one lattice distance right. Therefore, the transition shown in Figure 5 is from state [0,0,0,0] to state [-1,0,0,0]. Obviously, if all these states share the same value, they are equivalent to the ground state, but the entire molecule is translated laterally by the value; for example, [1,1,1,1] corresponds to the entire molecule moved by one lattice distance to the right. Figure 6 shows all the accessible basins for molecule A to move one lattice distance along the surface and the pathways and barriers separating them. Only the states for the molecule to move to the right are shown ([0,0,0,0][1,1,1,1]) because the basins for the molecule to move in the opposite direction ([0,0,0,0]-[-1,-1,-1,-1]) are simply a mirror image of these. The configuration of this molecule adsorbed in its alternative, rotated configuration is shown in Figure 4. Due to the different orientation of the molecule with respect to the surface rows, the mechanical properties this system and the commensurability of the molecule with the surface structure are different, resulting in a different potential energy surface. The accessible intermediate states are also mapped out and are shown in Table 1.

19580 J. Phys. Chem. C, Vol. 112, No. 49, 2008

Trevethan and Shluger

Figure 6. Diagram showing all the accessible states for molecule A to move one lattice distance to the right, with the values of the energy barriers between states shown.

TABLE 1: Energies of All Nonequivalent Transition Barriers (eV) transition

molecule A

molecule A (rot.)

molecule B

molecule C

(0,0,0,0 to -1,0,0,0) (-1,0,0,0 to 0,0,0,0) (-1,0,0,0 to -1,-1,0,0) (-1,-1,0,0 to -1,0,0,0) (-1,0,0,0 to -1,0,0,-1) (-1,0,0,-1 to -1,0,0,0)

0.53 0.16 0.40 0.29 0.60 0.17

0.38 0.16 0.45 0.35 0.50 0.18

0.64 0.38 0.36 0.64 0.59 0.27

0.45 0.18 0.61 0.38 0.23 0.40

We next consider the adsorption of molecule B, where the global minimum energy configuration has the same basic structure as that of molecule A, with the binding groups in the same positions along the rows, as shown in Figure 6. As with molecule A, the anthracene backbone straddles the oxygen row in approximately the same configuration, except that the backbone is lifted slightly farther from the surface becuase of the longer legs (the additional -CH2-). The same procedure for mapping out the potential energy basins and paths was then applied to this molecule. It was found that the same number and types of states are accessible, with the same structures (in terms of the binding group positions), but with the relative energies of these states and the separating barriers differing. For example, the minimum energy path to move from state [0,0,0,0] to state [-1,0,0,0] is shown in Figure 4, along with the same transition for molecule A. Here, it can be seen that the barrier to move from [0,0,0,0] to [-1,0,0,0] is larger (0.64 eV, as compared to 0.53 eV), and the barrier to move back from [-1,0,0,0] to [0,0,0,0] is also larger (0.38 eV, as compared to 0.16 eV). This can be explained by the fact that the longer legs of molecule B allow significantly more flexibility. The relative energies and barriers for molecule B, in comparison with molecule A, are listed in Table 1. For molecule C, which has a longer tetracene backbone, we again find that the minimum energy configuration has the same structural motif in terms of the binding groups, which is shown in Figure 8i. This molecule is considerably less commensurate with the surface structure than molecule A in the same orientation, as the distance between the binding groups attached to the molecule along the direction of the rows has been increased, which, as can be seen in the figure, distorts the structure of the tetracene backbone. This has the effect of pushing the binding groups out and away from the molecule

Figure 7. (i) Side view and (ii) end view of the ground-state configuration of molecule B.

and lowering the energy barrier to enter a leading leg state from the ground state to 0.45 eV. Figure 8ii shows molecule C in the [0,1,1,0] configuration, which is only 0.1 eV higher in energy than the ground state. The energy barriers for all the nonequivalent transitions for this system are listed in Table 1. Dynamical Evolution The thermally induced motion of the molecules on the surface will consist of the system’s moving between the potential energy basins described and determined in the previous section, socalled intramolecular transitions. The diffusion of the entire molecule, which is essentially the process of the molecule moving from the global minimum energy state to an equivalent state a single lattice separation away, will occur through a sequence of these intramolecular transitions. Since the barriers for each of the accessible intramolecular transitions are known, the real-time evolution of the system at a finite temperature can be evaluated using the kinetic Monte Carlo algorithm.16 In this algorithm, which is described in detail in ref 17, the rate for a

Diffusive Motion of Large Organic Molecules

J. Phys. Chem. C, Vol. 112, No. 49, 2008 19581

Figure 9. Snapshot of the real time evolution of the lateral position of molecule A.

Figure 8. (i) Side view of molecule C in its ground-state configuration. (ii) Side view of molecule C in the [0,1,1,0] configuration.

particular transition to occur is determined from the familiar Vineyard expression,18

r ) νe-β∆E

(1)

where ∆E is the energy barrier, ν is the attempt frequency prefactor, and β ) 1/kT. In the calculations performed here, we use the commonly applied value for the prefactor frequency of 1012 Hz. The actual frequencies for different transitions may differ from this and will depend on the dynamical coupling of the binding groups to both the surface and the rest of the molecule. However, in this initial model treatment where we are primarily concerned with understanding the different mechanisms of diffusion and the qualitative differences in the behavior of the molecules, this approximation is justified. The algorithm is implemented in a Fortran code, which stores all of the values for the energy barriers for a particular molecule-surface system in a large array, which is then used to determine the trajectory of the system on the potential energy surface using random numbers generated to determine both the time for a transition and the transition pathway. Each transition (and the time at which it occurs) is recorded as the simulation progresses. When the system moves into an equivalent ground state (a molecular transition), this is also recorded. The average time for the molecule to move between ground states is then used to determine the effective rate of diffusion of the entire molecule. All of the simulations presented in this paper have been performed using a temperature of 300 K unless otherwise stated. Molecule A We first examine the dynamics of molecule A, for which a snapshot of the evolution of the position along the surface as a function of time is shown in Figure 9 and an animation of the motion is included as Supporting Information. The average time for a single complete molecular transition in this case is 4.3 s, which leads to an effective diffusion rate of 0.12 Hz. This is a very slow rate of diffusion compared to the time scales of the intramolecular transitions and is due to the fact that on average,

44 000 intramolecular transitions are required for a single molecular translation to a neighboring ground state. This can be explained by referring to the potential energy surface of this molecule shown in Figure 6. Initially in the ground state, there are four equivalent states that the system can move to; that is, the movement of one of the four legs (the leading leg) to a neighboring basin (i.e., [-1,0,0,0], [0,1,0,0], [0,0,1,0], and [0,0,0,-1]). The barrier for each of these transitions is 0.53 eV, which will occur fairly quickly, on the order of a few milliseconds. Once the system has moved into one of these states, there are then three different states in each case that the molecule can move into: it can return to the same initial ground state with a barrier of 0.16 eV or move to one of the states depicted in the central column of Figure 6. The states depicted in the top and bottom of this column correspond to when the leg is in the same Ti row as the leading leg, follows the leading leg into the adjacent state (this state will be called the “walk” state, and there are four of them); the barrier for this transition is 0.40 eV. The state depicted in the center of this column occurs when the leg in the Ti row opposite the leading leg moves into a position alongside the leading leg (this state will be called the “crawl” state); the barrier for this transition is 0.60 eV. From these three different barriers, it is clear that the leading leg will follow the path back to the ground state in the vast majority of cases (in fact, with a probability of 0.9999). So when the system leaves its ground state and moves into one of the “leading leg” states, 9999 times out of 10 000, it will then return back to the initial state. One time out of the 10 000 it will enter the corresponding “walk” state, where it will then either return to the original “leading leg” state or move to a new “leading leg” state, both with barriers of 0.29 eV. If it moves into a new leading leg state, this will quickly fall (with a probability of 0.9999) into a new ground state, and the entire molecule will have moved one lattice distance. Therefore, on average, for every 20 000 transitions out of the ground state, only one will culminate in a complete molecular transition, resulting in the relatively high residence time of this molecule. The “crawl” state is visited 10 000 times less often that the “walk” state in this system because of the higher barrier (0.6 eV, as compared with 0.4 eV). When the temperature is increased, the “crawl” state is visited more frequently, but at 600 K, still only about 0.02 of successful molecular transitions will go through it. With molecule A adsorbed in its alternative configuration, with the backbone rotated 90° with respect to the rows of the surface, the dynamical evolution is significantly different as a result of the different potential energy surface, leading to an

19582 J. Phys. Chem. C, Vol. 112, No. 49, 2008 effective diffusion rate of 6 Hz, ∼50 times faster than in the original configuration. In this case, the barrier to return from the leading leg state is again much smaller than the barriers to move into the walk or crawl states; however, the barrier to leave the ground state is now only 0.38 eV, resulting in the much faster motion. In this configuration, the barriers to move into the walk and crawl states are quite similar, leading to 0.04 of the transitions occurring via the crawl state and 0.96 occurring via the walk state. The dependence on the rate of diffusion of molecule A on its orientation on the surface is an effect very similar to that observed in ref 19, where violet lander molecules were deposited on the anisotropic Cu(110) surface. By switching the orientation of the molecule with respect to the surface with a scanning tunnelling microscope tip, the diffusion rate is dramatically changed. The reason for this is again due to the different commensurability of the molecule to the surface in the different configurations. Molecule B We now examine the dynamical behavior of molecule B, adsorbed with the anthracene backbone parallel to the surface rows, which as described in the previous section has the same number and type of accessible states as molecule A. For this molecule, the average time of a complete molecular transition is 0.14 s, giving an effective diffusion rate of 3.5 Hz. This molecule then diffuses approximately 30 times faster than molecule A, even though all of the energy barriers for the intramolecular transitions are higher than in molecule A. This is due to the fact that now, on average, only 16 intramolecular transitions occur for a single molecular transition to a neighboring ground state. An animation of the real-time motion of molecule is included as Supporting Information. As before, the system in its ground state can move into one of four “leading leg” states, this time over a barrier of 0.64 eV, which has a transition rate 70 times lower than for the same process in molecule A. However, once in a leading leg state, the barrier to return to the ground state is now also much larger (see Figure 5) at 0.38 eV; larger, in fact, than the 0.36 eV barrier to move into the corresponding “walk” state. So now, even though the system leaves the ground state 70 times less often than molecule A, approximately one in every eight of these transitions will lead to the diffusion of the entire molecule. As before, the “crawl” state is negligibly occupied. Molecule C Moving to molecule C, this system again has the same number of states as before, but with a potential energy surface qualitatively different from molecule A because of the increased length of the aromatic backbone. For this molecule, the average time for a complete molecular transition is 0.000 14 s, giving an effective diffusion rate of 3.7 kHz, a rate ∼30 000 times faster than molecule A. The primary reason for this is that the energy barriers for many of the intramolecular transitions are significantly reduced as a result of the lower commensurability of the molecule with the positions of the Ti atoms in the surface. The distance between the groups along the molecule is 5.0 Å for the case of anthracene and 7.5 Å for the case of tetracene, whereas the distance between pairs of Ti atoms in the surface along the rows is 5.9 Å (see Figure 8). The barriers for the molecule to move from the ground state into one of the leading leg states is now only 0.45 eV, which is part of the reason this molecule is much more mobile. The other is that the barrier to then enter one of the “crawl” states (8ii) is

Trevethan and Shluger now only 0.23 eV, whereas the barrier to enter a “walk” state is 0.61 eV (this is also a result of the change in commensurability: the increased length of the molecule has made the “crawl” state much more favorable/accessible). Therefore, many more leading leg transitions will result in a complete molecular transition, for which, on average, 32 intramolecular transitions occur, all of which go through the “crawl” state. Discussion and Conclusion We have studied model systems consisting of three structurally different but chemically similar molecules adsorbed on an insulating surface. We investigated the real time evolution of these systems by first mapping out the accessible potential energy surface and then performing kinetic Monte Carlo simulations. This allowed us to determine the behavior of the system over long time scales in a computationally efficient way. The diffusion of these molecules occurs over time scales of milliseconds up to many seconds. More direct methods of studying diffusion, such as molecular dynamics, can only access time scales of up to nanoseconds. We have shown how these complex systems have a number of accessible and nonequivalent states and the relative rates of transition between these states that determines the mobility of the entire molecule. Very similar molecules can differ in their mobility on the surface by many orders of magnitude, which is due to the types of states that are accessible and the barriers that separate them. The potential energy landscape of this type of system is determined not only by the chemical interaction between the surface and the binding groups, but also by the structure of the entire molecule relative to the surface. With the molecules investigated in this study on the TiO2(110) surface, the differences in the potential energy surfaces were caused by the different lengths (and hence, flexibility) of the legs and overall length of the “backbone” structure. This was done both by changing the chemical structure of the molecules and by changing the orientation and adsorption structure of a molecule on the surface. From these results, it should be expected that limiting the freedom of movement and flexibility of the molecule should decrease the mobility, but at the same time, ensure that the binding parts of the molecule are commensurate with the surface. This is essentially the same design considerations reached for the diffusive properties of simpler two-legged molecules on the same surface.8 The systems investigated in this study are relatively simple and intuitive: any motion is restricted to one dimension, and the simple structural motif made it possible to easily understanding which states are likely to be accessible and identify them. This may not be the case, however, with other systems: a more isotropic surface would introduce many more possible states, as would more complex molecular structures that had less symmetry. In addition, it may not always be possible to intuitively identify all possible energy basins, and it may not be feasible to calculate all of the energy barriers, even with a classical atomistic model. In this case, more powerful methods for exploring the potential energy surface would be required, such a temperature-accelerated molecular dynamics20 or on-thefly KMC.21 In our calculations, we have assumed the pre-exponential frequency of the rate expression to be a commonly used value and the same for all the different transitions in each system. In this study, we believe that this approximation is sufficient to enable at least a qualitative comparison of the different molecular structures and to give insight into the diffusion mechanisms of this type of system. However, the actual

Diffusive Motion of Large Organic Molecules frequencies may differ and be dependent on both the local structure of the molecule and the dynamical coupling of the binding groups to the surface. A method such as temperatureaccelerated molecular dynamics20 would be required to determine the actual prefactors for each of the different processes. Another important limitation of this treatment is that the KMC method assumes the processes in the system are Markovian and that all the transitions are uncorrelated. If this is not the case, it may significantly change the behavior of the system. It may be the case that one “leg” has a significantly increased probability of being involved in a further transition, as opposed to any of the other legs, immediately after (within the correlation time) it has jumped, as a result of the local excitation. This effect will likely lead to higher residence times for the molecules investigated here (because of an increased chance of a leading leg state to return back to the ground state immediately after entering it). The effect on the intramolecular transitions of correlation effects will be affected by both the degree of coupling of the molecule to surface vibrational modes and how vibrational excitations are dissipated throughout the molecule. These effects could potentially be investigated through nonequilibrium molecular dynamics simulations. Acknowledgment. This work was supported through the EU FP6 project “PICO-inside”. The authors are grateful to M. Watkins for valuable discussions. Supporting Information Available: Animations (.mpg format) of the real time diffusion of molecules A and B. This information is available free of charge via the Internet at http:// pubs.acs.org.

J. Phys. Chem. C, Vol. 112, No. 49, 2008 19583 References and Notes (1) Dimitrakopoulos, C. D.; Purushothaman, S.; Kymissis, J.; Callegari, A.; Shaw, J. W. Science 1999, 283, 822. (2) Joachim, C.; Gimzewski, J. K.; Aviram, A. Nature 2000, 408, 541. (3) Meyer zu Heringdorf, F.-J.; Reuter, M. C.; Tromp, R. M. Nature 2001, 412, 517. (4) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press: Oxford, 2002. (5) Diebold, U. Surf. Sci. Rep. 2003, 48, 53. (6) Uetsuka, H.; Henderson, M. A.; Sasahara, A.; Onishi, H. J. Phys. Chem. B 2004, 108, 13706. (7) Gong, X. J. Phys. Chem. B 2006, 110, 2804. (8) Watkins, M.; Trevethan, T.; Sushko, M. L.; Shluger, A. L. J. Phys. Chem. C 2008, 112, 4226. (9) Trevethan, T.; Shluger, A. J. Phys. Chem. C 2007, 111, 15375. (10) Gundlach, L.; Szarko, J.; Socaciu-Siebert, L. D.; Neubauer, A.; Ernstorfer, R.; Willig, F. Phys. ReV. B 2007, 75, 125320. (11) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M., Jr.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1996, 118, 2309. (12) Bandura, A. V.; Kubicki, J. D. J. Phys. Chem. B 2003, 107, 11072. (13) Sushko, M. L.; Gal, A. Y.; Shluger, A. L. J. Phys. Chem. B 2006, 110, 4853. (14) Evarestov, R. A.; Bredow, Th.; Jug, K. Phys. Solid State 2001, 43, 1774. (15) Kantorovich, L.; Trevethan, T.; Foster, A. Self-Consistent Image Force Interaction code. http://www.cmmp.ucl.ac.uk/ lev/codes/SciFi/manual3-51/index.html. (16) Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. J. Comput. Phys. 1975, 17, 10. (17) Voter, A. F. In Introduction to the kinetic Monte Carlo method. Radiation Effects in Solids; Sickafus, K. F., Kotomin, E. A., Uberuaga, B. P., Eds.; Springer: Verlag: Berlin, 2007. (18) Vineyard, G. H. J. Phys. Chem. Solids 1957, 3, 121. (19) Otero, R.; Hummelink, F.; Sato, F.; Legoas, S. B.; Thostrup, P.; Laegsgaard, E.; Stensgaard, I.; Galvao, D. S.; Besenbacher, F. Nat. Mater. 2004, 3, 779. (20) Sorensen, M. R.; Voter, A. F. J. Chem. Phys. 2000, 112, 9599. (21) Henkelman, G.; Jonsson, H. J. Chem. Phys. 2001, 115, 9657.

JP806355M